JYU-Utils (empty) → 0.1.1
raw patch · 24 files changed
+1800/−0 lines, 24 filesdep +QuickCheckdep +arraydep +basesetup-changed
Dependencies added: QuickCheck, array, base, binary, bytestring, containers, deepseq, directory, filepath, haskell98, lazysmallcheck, mtl, mwc-random, parallel, process, random, stm, template-haskell, unix, zlib
Files
- JYU-Utils.cabal +29/−0
- LICENSE +19/−0
- Setup.lhs +4/−0
- Utils/BinaryInstances.hs +9/−0
- Utils/Darcs.hs +14/−0
- Utils/DynMap.hs +31/−0
- Utils/Exception.hs +24/−0
- Utils/File.hs +120/−0
- Utils/Function.hs +31/−0
- Utils/List.hs +264/−0
- Utils/Monad.hs +108/−0
- Utils/MonadRandom.hs +147/−0
- Utils/Numeric.hs +57/−0
- Utils/Parallel.hs +126/−0
- Utils/Point.hs +20/−0
- Utils/Rectangle.hs +120/−0
- Utils/Sampling.hs +36/−0
- Utils/SemanticEditors.hs +17/−0
- Utils/Shuffle.hs +68/−0
- Utils/ShuffleMWC.hs +75/−0
- Utils/Stream.hs +185/−0
- Utils/String.hs +57/−0
- Utils/Table.hs +88/−0
- Utils/Vector.hs +151/−0
+ JYU-Utils.cabal view
@@ -0,0 +1,29 @@+Name: JYU-Utils+Version: 0.1.1+Synopsis: Some utility functions for JYU projects+Description: List, function and monad utility functions. + Includes an old variant of MonadRandom that is used in + CV-package. ++ Pretty much the only reason to ever use this+ is to compile CV-package, and that is also just+ due legacy reasons. +Category: Utility+License: MIT+License-File: LICENSE+Author: Ville Tirronen+Maintainer: ville.tirronen@jyu.fi+Build-Type: Simple+Cabal-Version: >=1.4++Library+ Build-Depends: haskell98, base >= 3 && < 5, mtl >= 1.1.0, random >= 1.0.0, template-haskell, + process >= 1.0, QuickCheck >= 2.1, containers >= 0.2, stm >= 2.1, array >= 0.2, binary >= 0.5,+ zlib >= 0.5, bytestring >= 0.9, parallel >= 1.0, directory >= 1.0, filepath >= 1.1,+ unix >= 2.3, lazysmallcheck >= 0.5, mwc-random >= 0.8 && <0.9, deepseq >= 1.1 && < 2+ Exposed-modules: Utils.Darcs, Utils.Function, Utils.Numeric, Utils.Sampling, Utils.String, + Utils.DynMap,Utils.Parallel, Utils.SemanticEditors, + Utils.Table, Utils.File, Utils.List, Utils.Point, Utils.Shuffle, Utils.ShuffleMWC, Utils.Vector, + Utils.Monad, Utils.Rectangle, Utils.Stream, Utils.MonadRandom,+ Utils.BinaryInstances, Utils.Exception+
+ LICENSE view
@@ -0,0 +1,19 @@+Copyright (C) <year> by <copyright holders>++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in+all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN+THE SOFTWARE.
+ Setup.lhs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain
+ Utils/BinaryInstances.hs view
@@ -0,0 +1,9 @@+module Utils.BinaryInstances where++import Data.Binary++import Utils.Rectangle++instance (Binary a) => Binary (Rectangle a) where+ put (Rectangle a) = put a+ get = get >>= return.Rectangle
+ Utils/Darcs.hs view
@@ -0,0 +1,14 @@+{-# OPTIONS_GHC -XTemplateHaskell #-}+module Utils.Darcs where+import System.Process+import System.Exit+import Data.Maybe+import Language.Haskell.TH++getLastDarcsTag = do+ (exit,output,_) <- readProcessWithExitCode "darcs" ["show","tags"] [] + case exit of+ ExitSuccess -> return . Just . head . lines $ output+ _ -> return Nothing++thVersionTag = runIO getLastDarcsTag>>=litE . StringL . fromMaybe "Non versioned system"
+ Utils/DynMap.hs view
@@ -0,0 +1,31 @@+module Utils.DynMap where+import qualified Data.Map as Map+import Data.Dynamic+import Control.Monad.State+import Control.Monad+import Data.List(intercalate)++newtype PropertyMap = PM (Map.Map String (String,Dynamic))++emptyD :: PropertyMap+emptyD = PM Map.empty+insertD k v (PM dmap) = PM $ Map.insert k (show v,toDyn v) dmap+(PM dmap) #? k = fromDyn (snd $ dmap Map.! k) (error $ "TYPE ERROR ON "++k)+(PM dmap) #?? k = (fst $ dmap Map.! k) +(PM dmap) #/ k = (Map.delete k dmap) ++showD dmap = "{" ++ (intercalate ", " . map (\(k,(v,p)) -> k++": "++v) . Map.toList $ dmap) ++"}"++instance Show PropertyMap where+ show (PM pm) = "{"++ + (intercalate ", " . map (\(k,(v,_))-> k++":" ++v) $ Map.assocs pm)+ ++"}"++(=:) :: (Show a,Typeable a) => String -> a -> State PropertyMap ()+a =: b = do + x <- get+ put (insertD a b x)+ +properties st = snd $ runState st emptyD++
+ Utils/Exception.hs view
@@ -0,0 +1,24 @@+module Utils.Exception where+import Control.Exception as E+import System.Exit++inProcess a b = tagEM b a +tagEM op msg = op + `E.catch` + (\e -> error ("Error on "++msg+ ++" ("+ ++show (e::SomeException)+ ++")"))++tagE op msg = mapException tag op+ where + tag :: SomeException -> ErrorCall+ tag e = ErrorCall $ "Error in "++msg++" ("++show e++")"++onErrorTerminateWith op msg = op + `E.catch` + (\e -> putStrLn ("Error on "++msg+ ++" ("+ ++show (e::SomeException)+ ++")") + >> exitWith (ExitFailure 1))
+ Utils/File.hs view
@@ -0,0 +1,120 @@+{-# LANGUAGE FlexibleInstances, UndecidableInstances #-}+module Utils.File where+import System.Posix.Files(fileExist)+import System.Posix.Directory(readDirStream,openDirStream)+import System.FilePath.Posix+import System.Cmd+import System.Directory(createDirectoryIfMissing+ ,setCurrentDirectory+ ,getCurrentDirectory)+import Control.Exception as E(bracket_,catch,evaluate) +import qualified System.FilePath.Posix as Posix+import System.Directory+import System.IO++import Control.Parallel.Strategies+import Control.Exception+import IO+import qualified Data.ByteString.Char8 as B+import qualified Data.ByteString.Lazy as BS+import Codec.Compression.GZip++import Data.Binary++-- | Cache results of operation op under name 'fn'. +-- if correct file is found, its contents are used, otherwise 'op' is performed and it's results+-- are saved to the filed and returned from the call. ++class Cacheable a where+ writeCache :: FilePath -> a -> IO ()+ readCache :: FilePath -> IO a++box x = [x]++instance (Binary a) => Cacheable a+ where+ readCache fn = do+ (B.readFile fn >>= E.evaluate . decode . decompress . BS.fromChunks . box ) + `E.catch` (\err -> error $ "Error reading cache "++fn++": "+ ++show (err:: IOException))+ writeCache fn x = (BS.writeFile fn . compress . encode $ x)+ `E.catch` (\err -> error $ "Error writing cache "++fn++": "+ ++show (err:: IOException))++cached :: (Cacheable a) => FilePath -> IO a -> IO a+cached fn op = do+ let fn' = fn++".CACHE"+ e <- fileExist fn'+ if e+ then readCache fn'+ else do+ x <- op+ writeCache fn' x+ return x+++-- Check that file is not . or ..+notABackLink = (not . (flip elem) [".",".."] . takeFileName)++-- Get directory contents with path appended to them+getDirectoryContentsWithPath path = do+ dc <- getDirectoryContents path+ return $ map (path Posix.</>) dc++-- Return a filename that does not exists+genFileName base ext = genFileName' names+ where + genFileName' :: [FilePath] -> IO FilePath+ genFileName' (n:ns) = do + exists <- fileExist n+ if exists+ then genFileName' ns+ else return n + names :: [FilePath]+ names = [base ++ show x++ext | x <- [1..]]++-- return files from directory+getDirectoryList fp = openDirStream fp >>= getDSContent + where + getDSContent ds = do + x <- readDirStream ds+ if x == "" then return []+ else do + xs <- getDSContent ds+ return ((fp++x):xs)++-- Does file have extension `ext` +hasExt ext fp = all (uncurry (==)) $ zip (reverse ext) (reverse fp)++-- Retreive files from `path` that have extension `ext`+getFilesOfExt ext path = getDirectoryList path >>= return.filter (hasExt ext)++-- Perform operation `op` in directory `dir`. If directory does not exist+-- It will be created+inDirectory :: FilePath -> IO a -> IO a +inDirectory dir op = do+ createDirectoryIfMissing True dir+ currDir <- getCurrentDirectory+ Control.Exception.bracket_ + (setCurrentDirectory dir)+ (setCurrentDirectory currDir)+ op++-- Append text to `file` strictly and hax it so that parallel writes are forcefully done +strictPersistentAppendFile p file string = B.appendFile file (B.pack string)+ where + catch :: Int -> IOException -> IO ()+ catch 0 e = Control.Exception.throw e+ catch n e = strictPersistentAppendFile (n-1) file string++strictAppendFile file string = B.appendFile file (B.pack string)+--strictAppendFile file text = strictAppendFile' file text `demanding` rnf text+--strictAppendFile' outputfile text = Control.Exception.catch +-- (do +-- appendhandle <- openFile outputfile (AppendMode)+-- hPutStr appendhandle text+-- hFlush appendhandle+-- hClose appendhandle)+-- (\e ->fail "")++
+ Utils/Function.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE PatternGuards #-}+module Utils.Function where++both f (x,y) = (f x, f y)++with f = \x -> (x, f x)++under a b = \x y -> a (b x y)++-- Numerical functions++affine1d (toA,toB) (fromA,fromB) x = x1+ where+ x0 = (x - toA)/(toB-toA)+ x1 = x0*(fromB-fromA) + fromA++mkFst f a = (f a, a)+mkSnd f a = (a, f a)+++-- For Ord. Find a proper location for these++minBy :: (a -> a -> Ordering) -> a -> a -> a+minBy op a b | LT <- op a b = a+ | EQ <- op a b = a+ | GT <- op a b = b++maxBy :: (a -> a -> Ordering) -> a -> a -> a+maxBy op a b | LT <- op a b = b+ | EQ <- op a b = a+ | GT <- op a b = a
+ Utils/List.hs view
@@ -0,0 +1,264 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Utils.List where++import Data.List+import Data.Function+import Data.Maybe+import Control.Arrow ((&&&))+import qualified Data.Map as M+import Test.QuickCheck+++-- | Group list into indevidual pairs: [1,2,3,4] => [(1,2),(3,4)]. +-- Works only with even number of elements+pairs [] = []+pairs [x] = error "Non-even list for pair function"+pairs (x:y:xs) = (x,y):pairs xs++-- | Undo pairs function+fromPairs [] = []+fromPairs ((x,y):xs) = x:y:fromPairs xs++prop_pairsFromTo xs = even (length xs) ==> xs == fromPairs (pairs xs)++-- | Group list into pairs: [1,2,3] => [(1,2),(2,3)]. +-- Works with non null lists++pairs1 x = zip x (tail x)++-- | Undo pairs1 function+fromPairs1 [] = []+fromPairs1 [(x,y)] = [x,y]+fromPairs1 ((x,y):xs) = x:fromPairs1 (xs)++prop_pairsFromTo1 xs = length xs > 1 ==> xs == fromPairs1 (pairs1 xs)++crease op = map (uncurry op) . pairs1+creaseM op = sequence . (crease op)++ranks f xs = map fst $ rankBy f xs++rankBy f xs = map (\(rank,(orig,val)) -> (rank,val))+ . sortBy (compare`on`(fst.snd))+ . zip [1..] + . sortBy (f`on`snd) + . zip [1..] + $ xs++clusterBy :: Ord b => (a -> b) -> [a] -> [[a]]+clusterBy f = M.elems . M.map reverse . M.fromListWith (++)+ . map (f &&& return)+++groupItems b a items = map ( (b . head) &&& map a) + . groupBy ((==)`on` b)+ . sortBy (comparing b) $ items ++-- Assoc-list lookup with default value+lookupDef d a lst = fromMaybe d $ lookup a lst++-- get all consecutive pairs of a list: +--pairings "kissa"+-- => [('k','i'),('i','s'),('s','s'),('s','a')]++pairings [] = []+pairings [x,y] = [(x,y)]+pairings (x:y:ys) = (x,y):pairings (y:ys)++-- Perform an operation for each in lst+forEach fun lst = unfoldr op ([],lst)+ where+ op (start,[]) = Nothing+ op (start,a:as) = Just (start++(fun a):as+ ,(start++[a],as)) ++forPairs fun lst lst2 = map (map fst) + $ forEach (\(a,b)->(fun a b,b))+ $ zip lst lst2++-- +replicateList n l = concat $ replicate n l+--++concatZipNub (a:as) (b:bs) + | a == b = a:concatZipNub as bs+ | a /= b = a:b:concatZipNub as bs +concatZipNub [] _ = []+concatZipNub _ [] = []++histogram binWidth values = (map len grouped)+ where+ len x = (snap (head x), fromIntegral (length x)) + min = minimum values+ max = maximum values+ grouped = group sorted+ sorted = sort $ map snap values+ snap x = binWidth*(fromIntegral $ floor (x/binWidth))++binList binWidth op ivs = zip bins (map op values) + where+ values = map (map snd) grouped+ bins = map (fst.head) grouped + grouped = groupBy (\(a,_) (b,_) -> a == b ) sorted+ sorted = sortBy (comparing fst) $ map snapIndex ivs+ snapIndex (i,v) = (binWidth*(i`div`binWidth),v)+ +++-- Map numeric list so it becomes zeromean+zeroMean lst = map (\x -> x - mean) lst + where mean = average lst++-- Take n best elements according to fitnesses++takeNAccordingTo n (fitnesses,elements) = + take n+ $ sortBy (comparing fst)+ $ zip fitnesses elements++-- Zip two lists by selection function+select c = zipWith (\a b -> if c a b then a else b)++-- Take half++takeHalf lst = take (length lst `div` 2) lst++splitToNParts n lst | n <= 0 = error "splitToNParts n <= 0"+ | otherwise = takeLengths (lengths (length lst) n) lst+ where+ lengths len n = zipWith (+) (replicate n (len`div`n)) (replicate (len`mod`n) 1++repeat 0)++prop_splitEq n xs = n>0 ==> concat (splitToNParts n xs) == xs+prop_splitLen n xs = n>0 && n<= (length xs) ==> length (splitToNParts n xs) == n++-- Count elements that match predicate p+count p = foldl (\sum i -> if p i then sum+1 else sum) 0 ++-- Count frequencies of elements in list+frequencies lst = map (\x -> (head x,genericLength x)) $ group $ sort lst+normalizeFrequencies ls = map (\(a,b) -> (a,b/sum (map snd ls))) ls+-- Count average of list+average s = sum s / (genericLength $ s) ++-- Take n smallest given op+smallestBy op n lst = smallestBy' op n lst [] +smallestBy' op n [] o = o+smallestBy' op n (i:input) [] = smallestBy' op n input [i]+smallestBy' op n (i:input) output@(o:os) + = smallestBy' op n input (take n $ insertBy op i output)+-- (sloppily) Count median of list+median s | odd len = sorted !! middle+ | otherwise = ((sorted !! middle) + + (sorted !! (middle -1))) / 2+ where+ middle = len `div` 2+ sorted = sort s+ len = length s++takeTail n lst = reverse $ take n $ reverse lst++-- Count standard deviation of a list +stdDev l = sqrt (sum (map (\x -> (x - avg)^2) l) + / genericLength l)+ where avg = average l++-- Transform a list so that nth element is sum of n first elements+cumulate [] = []+cumulate values = tail $ scanl (+) 0 values++schwartzianTransform :: (Ord a,Ord b) => (a -> b) -> [a] -> [a]+schwartzianTransform f = map snd . sort . map (\x -> (f x, x))++sortVia f = map snd . sortBy cmp . map (\x -> (f x , x))+ where cmp (a1,a2) (b1,b2) = compare a1 b1++comparing p a b = compare (p a) (p b)++-- Pick element that has majority in the list+majority lst = head $ maximumBy (comparing length) $ group $ sort lst++-- Get all possible k-sized neighbourhoods in the list+getKNeighbourhoods k p = get (length p) pknot+ where + pknot = p++pknot+ get 0 p = []+ get i p = take k p:get (i-1) (tail p)++prop_headIdentical_KN n xs = 1 <= n && length xs >= 1 ==>+ map head (getKNeighbourhoods n xs)+ ==+ xs++-- Split a list to `l` length pieces.+splitToLength l lst = unfoldr split lst+ where+ split [] = Nothing+ split lst = Just (take l lst, drop l lst) ++-- Take n pieces of given lengths+--takeLengths :: [Int] -> [a] -> [[a]]+takeLengths [] lst = []+takeLengths (l:ls) lst = take l lst:takeLengths ls (drop l lst) ++prop_takeLen ls xs = all (>=0) ls && sum ls < length xs ==> length (takeLengths ls xs) == length ls+prop_takeLens ls xs = all (>=0) ls && sum ls < length xs ==> map length (takeLengths ls xs) == ls++-- From LicencedPreludeExts (hawiki)+splitBy :: (a->Bool) -> [a] -> [[a]]+splitBy _ [] = []+splitBy f list = first : splitBy f (dropWhile f rest)+ where+ (first, rest) = break f list++--splitBetween :: ((a,a) -> Bool) -> [a] -> [[a]]+splitBetween c acc [] = [reverse acc] +splitBetween c acc [a] = [reverse $ a:acc] +splitBetween c acc (a:b:cs) | c a b = (reverse $ a:acc):splitBetween c [] (b:cs)+ | otherwise = splitBetween c (a:acc) (b:cs)++-- split list into subsets matching predicate+tear op l = (filter (not.op) l, filter op l)++swapEverywhere a b = concat $ zipWith merge (inits a) (tails a)+ where+ merge i [] = []+ merge i (t:ts) = map (\x -> i++[x]++ts) b+++takeWhile2 op lst = reverse $ tw op [head lst] (tail lst)+ where+ tw _ l [] = []+ tw op l (x:xs) = if op (head l) x + then tw op (x:l) xs+ else l++applyMap val ops = map (\op -> op val) ops +applyMapM :: (Monad m) => a -> [a -> m b] -> m [b]+applyMapM val ops = mapM (\op -> op val) ops +changesM :: (Monad m) => [a -> m b] -> a -> m [b]+changesM = flip applyMapM++rollList (a:xs) = xs ++[a]+roll = rollList++mergeList a b = a ++ drop (length a) b++takeWhile1 test [] = []+takeWhile1 test (x:xs) | test x = x:takeWhile1 test xs+ | otherwise = [x]+++-- Modify each element in list with function that has knowledge of already+-- modified elements+editingMap f l = editingTrav f [] l++editingTrav fun [] l@(x:xs) = editingTrav fun [(fun l x)] xs+editingTrav fun a [] = reverse a+editingTrav fun ss l@(x:xs) = editingTrav fun+ (fun (reverse ss++l) x:ss)+ xs+++-- Rotations of list+rotate (x:xs) = xs++[x]+cycles x = take (length x) $ iterate rotate x
+ Utils/Monad.hs view
@@ -0,0 +1,108 @@+module Utils.Monad where+import Control.Monad+import Control.Parallel.Strategies( ($|) , rdeepseq)++-- Return result fully evaluated+bangReturn x = return $| rdeepseq $ x++-- Monadic if+conditional cond thenOp elseOp = do+ t <- cond+ if t then thenOp+ else elseOp++whenM cond thenOp = do+ t <- cond+ if t then thenOp+ else return ()+++pairM a b = do+ x <- a+ y <- b+ return (x,y)++sequenceWithPar init [] = init+sequenceWithPar init (op:ops) = sequenceWithPar (op init) ops++sequenceWithParM :: (Monad m) => a -> [a -> m a] -> m a+sequenceWithParM init [] = return init+sequenceWithParM init (op:ops) = do+ r <- op init+ sequenceWithParM r ops++scanJoinMaybe s ops = do+ sr <- s + cs sr ops [] + where + cs start (op:ops) res = do + r <- op start+ case r of+ Just r' -> cs r' ops (start:res)+ Nothing -> return $ reverse $ start:res+ cs start [] res = return $ reverse $ start:res ++scanJoin s ops = do+ sr <- s + cs sr ops [] + where + cs start (op:ops) res = do + r <- op start+ cs r ops (start:res)+ cs start [] res = return $ reverse $ start:res + ++-- chain ms = foldl1 (>=>) ms++doWhile :: (Monad m) => (a -> m (Maybe a)) -> a -> m a+doWhile op x = do+ r <- op x+ case r of + Just b -> do+ a <- doWhile op b+ return a+ Nothing -> return x++nTimes n op = sequence $ replicate n op+nTimes_ n op = sequence_ $ replicate n op++-- Keep performing action `op` until its result satisfies cond+untilM :: (Monad m) => (a -> Bool) -> (m a) -> m a+untilM cond op = do+ x <- op+ if cond x then return x else untilM cond op++repeatMLast n initial op = foldM join (initial) (replicate n op) >>= return.reverse+ where + join a b = b a+ ++repeatM n initial op = foldM join [initial] (replicate n op) >>= return.reverse+ where + join a b = do+ r <- b $ head a+ return $ r:a++iterateWhileM_ op cond x = recurse x + where+ recurse p = do+ r <- op p+ if cond r + then recurse r+ else return ()++iterateWhileM op cond x = recurse x [x] + where+ recurse p acc = do+ r <- op p+ if cond r + then recurse r (r:acc)+ else return $ r:acc++stepWhile op x = sw x [] + where+ sw p acc = do+ r <- op p+ case r of + Just b -> sw b (p:acc)+ Nothing -> p:acc
+ Utils/MonadRandom.hs view
@@ -0,0 +1,147 @@+{-# OPTIONS_GHC -fglasgow-exts #-}++-- Taken from haskell wiki: http://www.haskell.org/hawiki/MonadRandom++module Utils.MonadRandom (+ MonadRandom,+ getRandom,+ getRandomR,uniformRandomVector,uniformRandomVectorRS,+ gaussianVector,gaussianPerturbation,gaussianPerturbationR,gaussianRand,+ normRand,+ cauchyR,+ evalRandomT,+ runRandomT,+ evalRandomIO,+ evalRand,+ fromList,+ fromNonWeightedList,+ chooseAtRandom,+ randomList,+ randomRoll,+ withRandomPair,withRandomThree, withProbability,+ Rand, RandomT -- but not the data constructors+ ) where++import Random+import Control.Monad.State+import Control.Monad.Identity+import Data.List(genericLength)+import Utils.Monad+--import System.Random.Mersenne.Pure64++class (Monad m) => MonadRandom m where+ getRandom :: (Random a) => m a+ getRandomR :: (Random a) => (a,a) -> m a++newtype (RandomGen g) => RandomT g m a = RandomT { unRT :: StateT g m a }+ deriving (Functor, Monad, MonadTrans, MonadIO)++liftState :: (MonadState s m) => (s -> (a,s)) -> m a+liftState t = do v <- get+ let (x, v') = t v+ put v'+ return x++instance (Monad m, RandomGen g) => MonadRandom (RandomT g m) where+ getRandom = (RandomT . liftState) random+ getRandomR (x,y) = (RandomT . liftState) (randomR (x,y))++evalRandomT :: (RandomGen g, Monad m) => RandomT g m a -> g -> m a+evalRandomT x g = evalStateT (unRT x) g+runRandomT x g = runStateT (unRT x) g++evalRandomIO x = do+ g <- liftIO getStdGen+ (val,g') <- runStateT (unRT x) g+ liftIO $ setStdGen g'+ return val++--evalRandomMersenneIO seed x = do+-- let g = pureMT seed+-- (val,g') <- runStateT (unRT x) g+-- -- liftIO $ setStdGen g'+-- return val+--+++-- Boring random monad :)+newtype Rand g a = Rand { unRand :: RandomT g Identity a }+ deriving (Functor, Monad, MonadRandom)++evalRand :: (RandomGen g) => Rand g a -> g -> a+evalRand x g = runIdentity (evalRandomT (unRand x) g)++-- functions for performing common tasks.++uniformRandomVector l = replicateM l getRandom+uniformRandomVectorRS rngs = mapM getRandomR rngs++-- Generate normally distributed numbers using box-muller transformation.+-- Result is distribution with zero mean and deviation of 1.+-- Not numerically stable++gaussianRand :: (MonadRandom m, Random a, Floating a) => m [a]+gaussianRand = do+ x1 <- getRandomR (0,1)+ x2 <- getRandomR (0,1)+ return [sqrt (-2 * log x1 ) * cos (2*pi*x2)+ ,sqrt (-2 * log x1 ) * sin (2*pi*x2)]+ +gaussianVector sigma length = + replicateM (length `div` 2 + 1) gaussianRand+ >>= return.map (*sigma).take length.concat++gaussianPerturbationR sigmas lst = do+ mutation <- gaussianVector 1 + (length $ lst)+ return $ zipWith (+) lst $ zipWith (*) sigmas mutation++gaussianPerturbation sigma lst = do+ mutation <- gaussianVector sigma + (length $ lst)+ return $ zipWith (+) lst mutation++normRand mu sigma = gaussianRand >>= return.(\x -> mu + sigma*x).head++-- fromList :: (MonadRandom m) => [(a,Rational)] -> m a+fromList [] = error "MonadRandom.fromList called with empty list"+fromList [(x,_)] = return x+fromList xs = do + let s = fromRational $ sum (map snd xs) -- total weight+ cs = scanl1 (\(x,q) (y,s) -> (y, s+q)) xs + -- cumulative weight+ p <- liftM toRational $ getRandomR (0.0,s::Double)+ return $ fst $ head $ dropWhile (\(x,q) -> q < p) cs++fromNonWeightedList l = do+ index <- getRandomR (0,genericLength l -1)+ return $ l !! index++chooseAtRandom (a,b) = do+ x <- getRandomR (0,1:: Int)+ if x >0 then return a else return b++randomList items n = replicateM n $ fromNonWeightedList items+randomRoll lst = do+ i <- getRandomR (0,length lst-1)+ return $ drop i lst ++ take i lst++withRandomPair op population = do+ a <- fromNonWeightedList population + b <- fromNonWeightedList population + op a b++withRandomThree op population = do+ a <- fromNonWeightedList population + b <- fromNonWeightedList population + c <- fromNonWeightedList population + op a b c++withProbability :: MonadRandom m => Double -> m a -> m a -> m a+withProbability p op op2 = do+ u <- getRandom+ if u<p then op else op2++cauchyR x0 gamma = do+ y <- getRandomR (0,1::Double)+ return $ gamma*tan(pi*(y-0.5))+x0
+ Utils/Numeric.hs view
@@ -0,0 +1,57 @@+module Utils.Numeric where+--import Numeric.LinearAlgebra +import Control.Monad+import Utils.MonadRandom+import Utils.Shuffle++-- Utilities for handling of numerical lists+inBounds bounds vector = all (\((a,b),x) -> a < x && x < b) $ zip bounds vector+saturateRealsFit f bounds (_,x) = (f (zipWith saturate bounds x)+ ,zipWith saturate bounds x)++saturateReals bounds x = zipWith saturate bounds x+saturateRealsM bounds x = return $ saturateReals bounds x++saturate (lb,ub) x = min (max x lb) ub++widths :: Num a => [(a,a)] -> [a]+widths = map (\(a,b) -> abs (a-b)) +++clean x | isNaN x || isInfinite x = 0+ | otherwise = x++{-scaleToUnitRows = fromRows . map (\v -> (1/sqrt (v<.>v)) .* v ) . toRows++gramSchmidt = fromRows . gs [] . toRows+gs :: [Vector Double] -> [Vector Double] -> [Vector Double]+gs [] (v:vs) = gs [v] vs+gs us (v:vs) = gs (vr:us) vs+ where vr = foldl1 (-) $ v:[proj ui v | ui <- us]+gs us [] = reverse us -- the silly bit+proj u v = ((v <.> u) / (u <.> u)) .* u+-}++-- Generate truly orthogonal random projection matrix. +{-+mkRandomProjMatrix1 :: (MonadRandom m) => Int -> Int -> m (Matrix Double)+mkRandomProjMatrix1 fromD toD = do+ numbers <- replicateM (fromD*toD) $ getRandomR (-1,1::Double)+ shuffled <- doShuffle . map (\v -> (1/sqrt (v<.>v)) .* v ) + . gs [] . toRows $ (fromD><toD) numbers+ return $ fromRows shuffled++-- Generate gaussian matrix+mkRandomProjMatrix2 :: (MonadRandom m) => Int -> Int -> m (Matrix Double)+mkRandomProjMatrix2 fromD toD = do+ numbers <- gaussianVector (1::Double) (fromD*toD) + return $ (fromD><toD) numbers++-- Generate RP matrix according to Achlioptas+mkRandomProjMatrix3 :: (MonadRandom m) => Int -> Int -> m (Matrix Double)+mkRandomProjMatrix3 fromD toD = do+ numbers <- replicateM (fromD*toD) $ MonadRandom.fromList [(-1,1/6)+ ,(0 ,2/3)+ ,(1 ,1/6)]+ return $ (fromD><toD) numbers+-}
+ Utils/Parallel.hs view
@@ -0,0 +1,126 @@+module Utils.Parallel where+import Data.IORef+import Control.Concurrent+import Control.Concurrent.STM+import Control.Exception as E+import Control.Monad++import Utils.Monad+import Utils.File++-- Module for running simple tasks in parallel+-- Fork with indicator when the action completes+-- myForkOS :: IO () -> IO (MVar ())+myForkIO io = do+ mvar <- newEmptyMVar + forkIO (io `finally` (putMVar mvar ()) )+ return mvar++-- Run IO () actions in parallel. It is assumed that the return values are+-- stored elsewhere++tvarTake x n = do + l <- readTVar x+ let (res,rem) = (take n l,drop n l)+ writeTVar x rem+ return res++tvarPop x = do + l <- readTVar x+ writeTVar x $ tail l+ return $ head l++tvarPush x t = do + l <- readTVar x+ writeTVar x $ l ++ [t]+ return $ head l++liftT f = \var -> do+ v <- readTVar var+ f v++data TaskRunnerStrategy = Persistent | Failing deriving(Eq,Show) ++taskRunner :: TaskRunnerStrategy -> TVar Bool -> TVar [IO ()] -> IO ()+taskRunner s terminated tasks = complete (Just $ return ())+ where+ catch :: IO () -> IOException -> IO ()+ catch t e = case s of + Persistent -> do+ atomically (tvarPush tasks t)+ strictAppendFile "PAR-RUN-ERRORS" ("persisting: "++show e++"\n")+ `E.catch` (\err -> error $ "Error adding to PAR-RUN-ERRORS "+ ++show (err:: IOException)++". Tried to add "+ ++show e)+ + Failing -> strictAppendFile "PAR-RUN-ERRORS" ("failing: "++show e++"\n") + `E.catch` (\err -> error $ "Error adding to PAR-RUN-ERRORS "+ ++show (err:: IOException)++". Tried to add "+ ++show e)++ + complete (Just task) = do + task `E.catch` (catch task)+ newTask <- next+ complete newTask+ complete Nothing = atomically $ writeTVar terminated True++ next = atomically $ do+ anyLeft <- (liftT (return.not.null)) tasks+ if anyLeft then tvarPop tasks + >>= \x -> return $ Just x+ else return Nothing++parRun :: TaskRunnerStrategy -> Int -> [IO ()] -> IO ()+parRun s poolSize tasks = do+ t <- atomically $ newTVar tasks+ terminators <- replicateM poolSize $ do+ terminator <- atomically $ newTVar False+ myForkIO (taskRunner s terminator t)+ return terminator+ atomically $ do+ terminated <- mapM readTVar terminators >>= return.and -- (liftT (return.null)) t+ check terminated+++parRunWithMonitor :: TaskRunnerStrategy -> Int -> (TVar [IO ()] -> [TVar Bool] -> IO ()) -> [IO ()] -> IO ()++parRunWithMonitor s poolSize monitor tasks = do+ t <- atomically $ newTVar tasks+ terminators <- replicateM poolSize $ do+ terminator <- atomically $ newTVar False+ myForkIO (taskRunner s terminator t)+ return terminator+ forkIO $ (monitor t terminators)+ atomically $ do+ terminated <- mapM readTVar terminators >>= return.and -- (liftT (return.null)) t+ check terminated++--startPool :: Int -> [IO ()] -> IO ()+--startPool poolSize tasks = do+-- t <- atomically $ newTVar tasks+-- pool <- replicateM poolSize (myForkIO (tasks t))+-- atomically $ do+-- terminated <- (liftT (return.null)) t+-- check terminated++-- Simple Concurrent mapping+++fork1 :: (a -> IO b) -> a -> IO (MVar b)+fork1 f x =+ do+ cell <- newEmptyMVar+ forkIO (do { result <- f x; putMVar cell result })+ return cell++fork :: (a -> IO b) -> [a] -> IO [MVar b]+fork f = mapM (fork1 f)++joinMVars :: [MVar b] -> IO [b]+joinMVars = mapM takeMVar++forkAndJoin :: (a -> IO b) -> [a] -> IO [b]+forkAndJoin f xs = (fork f xs) >>= joinMVars++
+ Utils/Point.hs view
@@ -0,0 +1,20 @@+{-# LANGUAGE TypeSynonymInstances #-}+module Utils.Point where++type Pt a = (a,a)++instance Num a => Num (Pt a) where+ (x1,x2) + (y1,y2) = (x1+y1,x2+y2) + (x1,x2) * (y1,y2) = (x1*y1,x2*y2) + (x1,x2) - (y1,y2) = (x1-y1,x2-y2) + negate (x1,x2) = (negate x1,negate x2) + abs (x1,x2) = (abs x1, abs x2) -- Not really mathematical, but the type must be same+ signum (x1,x2) = (signum x1, signum x2) -- Ditto+ fromInteger x = (fromInteger x,fromInteger x) -- As well++norm2 :: (Num a) => Pt a -> a+norm2 (a,b) = a*a+b*b++norm = sqrt . norm2++(a,b) >/ (c,d) = (a `div` c, b `div` d)
+ Utils/Rectangle.hs view
@@ -0,0 +1,120 @@+module Utils.Rectangle where+import Test.LazySmallCheck++import Utils.Point+import Control.DeepSeq++newtype Rectangle a = Rectangle ((a,a),(a,a)) deriving (Eq,Show)++a `s` b = rnf a `seq` b ++instance (NFData a) => NFData (Rectangle a) where+ rnf (Rectangle ((a,b),(c,d))) = (a `s` b `s` c `s` d) `seq` ()++left (Rectangle ((x,y),(w,h))) = x+right (Rectangle ((x,y),(w,h))) = x+w+top (Rectangle ((x,y),(w,h))) = y+bottom (Rectangle ((x,y),(w,h))) = y+h+topLeft (Rectangle ((x,y),(w,h))) = (x,y)+topRight (Rectangle ((x,y),(w,h))) = (x+w,y) +bottomLeft (Rectangle ((x,y),(w,h))) = (x,y+h) +bottomRight (Rectangle ((x,y),(w,h))) = (x+w,y+h) +vertices r = [topLeft r, topRight r, bottomLeft r, bottomRight r]+rSize (Rectangle ((x,y),(w,h))) = (w,h)+rArea r = let (w,h) = rSize r in (w*h)++-- TODO: Add documentation #Cleanup++instance (Num a, Ord a , Serial a) => Serial (Rectangle a) where+ series = cons4 $ \a b c d -> mkRectangle (a,b) (c,d)++-- | Create rectangle around point (x,y)+around (x,y) (w,h) = mkRectangle (x', y') (w,h)+ where (x',y') = (x-(w/2),y-(h/2))++mkRectangle (x,y) (w,h) = Rectangle ((x-negW,y-negH),(abs w,abs h))+ where+ negH | h<0 = abs h+ | h>=0 = 0+ negW | w<0 = abs w+ | w>=0 = 0++mkRectCorners (x1,y1) (x2,y2) = Rectangle ((x,y),(w,h))+ where+ x = min x1 x2+ y = min y1 y2+ w = abs (x1-x2)+ h = abs (y1-y2)++prop_Corners :: (Int,Int) -> (Int,Int) -> Bool+prop_Corners p w = mkRectCorners p (p+w) == mkRectangle p w++mkRec = uncurry mkRectangle++-- | Return rectangle r2 in coordinate system defined by r1+inCoords r1 r2@(Rectangle (pos,size)) = Rectangle (pos-topLeft r1,size )++-- | Return a point in coordinates of given rectangle+inCoords' r1 pt = pt - topLeft r1++-- | Adjust the size of the rectangle to be divisible by 2^n.+enlargeToNthPower n (Rectangle ((x,y),(w,h))) = Rectangle ((x,y),(w2,h2))+ where+ (w2,h2) = (pad w, pad h)+ pad x = x + (np - x `mod` np)+ np = 2^n++intersection r1 r2 + = mkRectCorners (max (left r1) (left r2)+ ,max (top r1) (top r2))+ (min (right r1) (right r2)+ ,min (bottom r1) (bottom r2))++propIntersectionArea r1 r2 + = (intersects r1 r2) + ==> rArea (intersection r1 r2) <= rArea r1 &&+ rArea (intersection r1 r2) <= rArea r2++propIntersectionCommutes r1 r2 + = (intersects r1 r2) + ==> (intersection r1 r2) == (intersection r2 r1) ++intersects rect1 rect2 + = intersect1D (left rect1, right rect1) (left rect2, right rect2) && + intersect1D (top rect1, bottom rect1) (top rect2, bottom rect2)++contains a b = left a <= left b + && top a <= top b + && bottom a >= bottom b+ && right a >= right b++intersect1D (x,y) (u,w) = + not $ (x < min u w && y < min u w) || (x > max u w && y > max u w) ++prop_intersect1DCommutes a b + = intersect1D a b == intersect1D b a++prop_intersectsCommutes sa@(_,(s1,s2)) sb@(b,(s3,s4)) + = intersects (mkRec sa) (mkRec sb) == intersects (mkRec sb) (mkRec sa)++-- | Create a tiling of a rectangles. +tile tilesize overlap r = [mkRectangle ((x,y)-overlap) tilesize + | x <- [startx,startx+fst tilesize..endx]+ , y <- [starty,starty+fst tilesize..endy] ]+ where+ startx = left r-fst overlap+ starty = top r-snd overlap+ endx = right r+fst overlap+ endy = bottom r+snd overlap++-- | Scale a rectangle+scale (a,b) (Rectangle ((x,y),(s1,s2))) + = mkRectangle (round (a*fromIntegral x),round (b*fromIntegral y))+ (round (a*fromIntegral s1),round (b*fromIntegral s2))+++toInt (Rectangle (p, s)) + = Rectangle (both round p + ,both round s)+ where both f (a,b) = (f a , f b)+
+ Utils/Sampling.hs view
@@ -0,0 +1,36 @@+module Utils.Sampling where++import Utils.SemanticEditors+import Control.Arrow+import Control.Monad+import Data.List+import Utils.List+import Utils.MonadRandom+import Utils.Shuffle++-- This module is for separating sampling from MonadRandom in hopes of some day replacing it with+-- the standard version.++-- Uniform probability selection+selectRandom n l = replicateM n $ fromNonWeightedList l++-- Stochastic universal sampling+-- Results are randomized order+doSUS no weightedList = do+ x <- getRandomR (0,1) -- random spin+ doShuffle (sus weightedList no x)++sus :: [(Double,a)] -> Int -> Double -> [a]+sus elements no offsetRatio = sus' pts es [] + where+ offset = offsetRatio*armSpan+ armSpan = totalWeight/fromIntegral no+ es = ((inzip first (scanl1 (+)).cycle.sortBy (comparing (negate.fst))) elements)+ pts = [min totalWeight (offset+i) | i <- take no [0,armSpan..]]+ totalWeight = sum.map fst $ elements++sus' [] _ acc = acc+sus' a@(p:pts) r@((w,e):es) acc + | p <= w = sus' pts r (e:acc)+ | p > w = sus' a es acc+ | otherwise = sus' pts r acc
+ Utils/SemanticEditors.hs view
@@ -0,0 +1,17 @@+module Utils.SemanticEditors where+import Data.Array+import qualified Data.Map as Map+import Control.Arrow+-- Semantic editors, maybe to move somewhere else++ -- Array element+element i = \e arr -> arr // [(i,e (arr ! i))]+ +-- Finite Map +value k = \e map -> Map.adjust e k map++-- Lists of pairs -- Not a combinator+inzipmap sel fs = uncurry zip . sel (fmap fs) . unzip+inzip sel fs = uncurry zip . sel (fs) . unzip++
+ Utils/Shuffle.hs view
@@ -0,0 +1,68 @@+{-#LANGUAGE BangPatterns#-}+module Utils.Shuffle(shuffle,doShuffle) where+import Utils.MonadRandom++-- A complete binary tree, of leaves and internal nodes.+-- Internal node: Node card l r+-- where card is the number of leaves under the node.+-- Invariant: card >=2. All internal tree nodes are always full.+data Tree a = Leaf a | Node Int (Tree a) (Tree a) deriving Show++fix f = g where g = f g -- The fixed point combinator++-- Convert a sequence (e1...en) to a complete binary tree++build_tree = (fix grow_level) . (map Leaf)+ where+ grow_level self [node] = node+ grow_level self l = self $ inner l+ + inner [] = []+ inner [!e] = [e]+ inner (e1:e2:rest) = (join e1 e2) : inner rest+ + join l@(Leaf _) r@(Leaf _) = Node 2 l r+ join l@(Node ct _ _) r@(Leaf _) = Node (ct+1) l r+ join l@(Leaf _) r@(Node ct _ _) = Node (ct+1) l r+ join l@(Node ctl _ _) r@(Node ctr _ _) = Node (ctl+ctr) l r+++-- given a sequence (e1,...en) to shuffle, and a sequence+-- (r1,...r[n-1]) of numbers such that r[i] is an independent sample+-- from a uniform random distribution [0..n-i], compute the+-- corresponding permutation of the input sequence.++doShuffle elements = do+ is <- sequence $ [getRandomR (0,length elements-i-1) | i<-[0..length elements-2]]+ return $ (shuffle elements is)+++shuffle elements rseq = shuffle1' (build_tree elements) rseq+ where+ shuffle1' (Leaf e) [] = [e]+ shuffle1' tree (r:r_others) = + let (b,rest) = extract_tree r tree+ in b:(shuffle1' rest r_others)+ + -- extract_tree n tree+ -- extracts the n-th element from the tree and returns+ -- that element, paired with a tree with the element+ -- deleted.+ -- The function maintains the invariant of the completeness+ -- of the tree: all internal nodes are always full.+ -- The collection of patterns below is deliberately not complete.+ -- All the missing cases may not occur (and if they do,+ -- that's an error.+ extract_tree 0 (Node _ (Leaf e) r) = (e,r)+ extract_tree 1 (Node 2 (Leaf l) (Leaf r)) = (r,Leaf l)+ extract_tree !n (Node c (Leaf l) r) =+ let (e,new_r) = extract_tree (n-1) r+ in (e,Node (c-1) (Leaf l) new_r)+ extract_tree n (Node n1 l (Leaf e)) + | n+1 == n1 = (e,l)+ + extract_tree n (Node c l@(Node cl _ _) r) + | n < cl = let (e,new_l) = extract_tree n l+ in (e,Node (c-1) new_l r)+ | otherwise = let (e,new_r) = extract_tree (n-cl) r+ in (e,Node (c-1) l new_r)
+ Utils/ShuffleMWC.hs view
@@ -0,0 +1,75 @@+{-#LANGUAGE BangPatterns, ScopedTypeVariables#-}+module Utils.ShuffleMWC(shuffle,doShuffle) where+import System.Random.MWC+import Control.Monad(replicateM)+import Control.Applicative((<$>))+import Data.List(genericLength)++-- A complete binary tree, of leaves and internal nodes.+-- Internal node: Node card l r+-- where card is the number of leaves under the node.+-- Invariant: card >=2. All internal tree nodes are always full.+data Tree a = Leaf a | Node Int (Tree a) (Tree a) deriving Show++fix f = g where g = f g -- The fixed point combinator++-- Convert a sequence (e1...en) to a complete binary tree++build_tree = (fix grow_level) . (map Leaf)+ where+ grow_level self [node] = node+ grow_level self l = self $ inner l+ + inner [] = []+ inner [!e] = [e]+ inner (e1:e2:rest) = (join e1 e2) : inner rest+ + join l@(Leaf _) r@(Leaf _) = Node 2 l r+ join l@(Node ct _ _) r@(Leaf _) = Node (ct+1) l r+ join l@(Leaf _) r@(Node ct _ _) = Node (ct+1) l r+ join l@(Node ctl _ _) r@(Node ctr _ _) = Node (ctl+ctr) l r+++inRange gen (a,b) = do+ i :: Float <- uniform gen+ return (a+i*(b-a))++doShuffle gen elements = do+ is <- sequence $ [floor <$> inRange gen (0,fromIntegral $ length elements-i)+ | i<-[0..length elements-2]]+ return $ (shuffle elements is)++-- given a sequence (e1,...en) to shuffle, and a sequence+-- (r1,...r[n-1]) of numbers such that r[i] is an independent sample+-- from a uniform random distribution [0..n-i], compute the+-- corresponding permutation of the input sequence.++shuffle elements rseq = shuffle1' (build_tree elements) rseq+ where+ shuffle1' (Leaf e) [] = [e]+ shuffle1' tree (r:r_others) = + let (b,rest) = extract_tree r tree+ in b:(shuffle1' rest r_others)+ + -- extract_tree n tree+ -- extracts the n-th element from the tree and returns+ -- that element, paired with a tree with the element+ -- deleted.+ -- The function maintains the invariant of the completeness+ -- of the tree: all internal nodes are always full.+ -- The collection of patterns below is deliberately not complete.+ -- All the missing cases may not occur (and if they do,+ -- that's an error.+ extract_tree 0 (Node _ (Leaf e) r) = (e,r)+ extract_tree 1 (Node 2 (Leaf l) (Leaf r)) = (r,Leaf l)+ extract_tree !n (Node c (Leaf l) r) =+ let (e,new_r) = extract_tree (n-1) r+ in (e,Node (c-1) (Leaf l) new_r)+ extract_tree n (Node n1 l (Leaf e)) + | n+1 == n1 = (e,l)+ + extract_tree n (Node c l@(Node cl _ _) r) + | n < cl = let (e,new_l) = extract_tree n l+ in (e,Node (c-1) new_l r)+ | otherwise = let (e,new_r) = extract_tree (n-cl) r+ in (e,Node (c-1) l new_r)
+ Utils/Stream.hs view
@@ -0,0 +1,185 @@+module Utils.Stream where+import Control.Applicative+import Control.Monad+import Control.Monad.Trans++-- This module provides Monadic streams.++-- | Stream of monadic values+data Stream m a = Terminated | Value (m (a,Stream m a))++-- | Attaching side effects+sideEffect :: (Monad m) => (a -> m ()) -> Stream m a -> Stream m a+sideEffect p Terminated = Terminated+sideEffect p (Value next) = Value renext+ where + renext = do+ (r,n) <- next+ p r+ return (r,sideEffect p n)++-- | Repeating stream+listToStream [] = Terminated+listToStream (l:lst) = Value (return (l,listToStream lst))+repeatS x = Value (return (x,repeatS x))+repeatSM x = sequenceS (repeatS x) +-- | Create a stream by iterating a monadic action+iterateS op n = Value cont+ where + cont = do+ r <- op n+ return $ (n,iterateS op r) ++-- | Pure and monadic left fold over a stream+foldS op i Terminated = return i+foldS op i (Value xs) = xs >>= \(x,xn) -> foldS op (op i x) xn++foldSM op i Terminated = return i+foldSM op i (Value xs) = xs >>= \(x,xn) -> op i x >>= \opix -> foldSM op opix xn++-- | Merge two (time)streams+time (a,_) = a+value (_,a) = a+mergeTimeStreams starta startb a b = mergeE (starta,startb) (mergeS a b)+mergeTimeStreamsWith sa sb op a b = fmap (\(t,(a,b)) -> (t,(op a b))) $ mergeTimeStreams sa sb a b+mergeManyW starts op streams = snd $ foldl1 (\(s,m) (s1,n) -> ((op s s1),mergeTimeStreamsWith s s1 op m n)) (zip starts streams)++mergeS Terminated _ = Terminated+mergeS _ Terminated = Terminated+mergeS _ Terminated = Terminated+mergeS (Value xs) (Value ys) = Value renext+ where+ renext = do+ (x,xn) <- xs+ (y,yn) <- ys+ case compare (time x) (time y) of+ LT -> return (L x,mergeS xn (push y yn))+ EQ -> return (B (time x,(value x,value y)),mergeS xn yn)+ GT -> return (R y,mergeS (push x xn) yn)++data LRB a b c = L a | B b | R c deriving (Show)++mergeE _ Terminated = Terminated+mergeE (l,r) (Value xs) = Value renext+ where+ renext = do+ (x,xn) <- xs+ case x of+ L (t,a) -> return ((t,(a,r)),mergeE (a,r) xn)+ B (t,(a,b)) -> return ((t,(a,b)),mergeE (a,b) xn)+ R (t,b) -> return ((t,(l,b)),mergeE (l,b) xn)++push x Terminated = Value (return (x,Terminated))+push x xs = Value (return (x,xs))+ ++-- | Map over a stream+instance (Monad m) => Functor (Stream m) where+ fmap _ Terminated = Terminated+ fmap f (Value next) = Value renext+ where + renext = do+ (r,n) <- next+ return (f r,fmap f n)++instance (Monad m) => Applicative (Stream m) where+ pure f = repeatS f+ Terminated <*> _ = Terminated+ _ <*> Terminated = Terminated+ (Value a) <*> (Value b) = Value renext + where + renext = do+ (fun,anext) <- a+ (br,bnext) <- b+ return (fun br,anext<*>bnext)+ +zipS a b = (,) <$> a <*> b+ +--+sequenceS :: (Monad m) => Stream m (m a) -> (Stream m a)+sequenceS Terminated = Terminated+sequenceS (Value next) = Value $ do+ (op,n) <- next+ r <- op + return (r,sequenceS n)++mapMS :: (Monad m) => (a -> m b) -> Stream m a -> Stream m b+mapMS op s = sequenceS . fmap op $ s++-- |Drop elements from the stream. Due to stream structure, this operation cannot+-- fail gracefully when dropping more elements than what is found in the stream+dropS :: (Monad m) => Int -> Stream m a -> Stream m a+dropS _ Terminated = Terminated+dropS n next = Value renext+ where+ drop 0 s = return s+ drop _ Terminated = return Terminated+ drop n (Value next) = do+ (r,ne) <- next + drop (n-1) ne+ renext = do+ r <- drop n next+ case r of+ Terminated -> error "Not enough elements to drop"+ Value x -> x++takeS :: (Monad m) => Int -> Stream m a -> Stream m a+takeS _ Terminated = Terminated+takeS n (Value next) = Value renext+ where+ renext = do+ (r,ne) <- next+ if n<1 then return (r,Terminated)+ else return (r,takeS (n-1) ne)++takeWhileS _ Terminated = Terminated+takeWhileS c (Value next) = Value renext+ where+ renext = do+ (r,ne) <- next+ if not . c $ r then return (r,Terminated)+ else return (r,takeWhileS c ne)++consS a Terminated = Value (return (a, Terminated))+consS a s = Value (return (a, s))++-- pairS is safe only for infinite streams+pairS :: (Monad m) => Stream m a -> Stream m (a,a)+pairS Terminated = Terminated+pairS (Value next) = Value renext+ where+ renext = do+ (val1,nexts2) <- next+ case nexts2 of+ Terminated -> return (undefined,Terminated)+ (Value next2) -> do (val2,next3) <- next2+ return ((val1,val2),pairS (consS val2 next3))+++terminateOn :: (Monad m) => (a -> Bool) -> Stream m a -> Stream m a+terminateOn cond Terminated = Terminated+terminateOn cond (Value next) = Value renext+ where+ renext = do+ (r,n) <- next+ if cond r then return (r,Terminated)+ else return (r,terminateOn cond n)++runStream Terminated = return []+runStream (Value s) = do+ (n,next) <- s+ r<-runStream next+ return (n:r)++runStream_ Terminated = return ()+runStream_ (Value s) = do+ (n,next) <- s+ runStream_ next++runLast l Terminated = return l+runLast l (Value s) = do+ (n,next) <- s+ runLast n next++runLast1 s = runLast (error "Empty Stream") s+
+ Utils/String.hs view
@@ -0,0 +1,57 @@+{-# OPTIONS -fglasgow-exts #-}+module Utils.String where++import Numeric+import Test.QuickCheck++-- Utilities for printing numbers+ +showPercentage n = showFFloat (Just 2) n ""+++-- TODO: Unify the three below+columns :: [(String,String)] -> String+columns pairs = unlines $ cd+ where+ fstLen = maximum $ map (length.fst) pairs + cd = map (\(a,b) -> padToR ' ' fstLen a++"\t"++b) pairs++columnS :: (Show a) => [(String,a)] -> String+columnS pairs = unlines $ cd+ where+ fstLen = maximum $ map (length.fst) pairs + cd = map (\(a,b) -> padToR ' ' fstLen a++"\t"++show b) pairs++columnsBy show pairs = unlines $ cd+ where+ fstLen = maximum $ map (length.fst) pairs + cd = map (\(a,b) -> padToR ' ' fstLen a++"\t"++show b) pairs++-- Pad string/list `s` to length `width` using `padding`+padToR pad width s = s ++ padding + where + padding = replicate n pad+ n = max 0 (width - length s)+padTo pad width s = padding ++ s + where + padding = replicate n pad+ n = max 0 (width - length s)++prop_padLength (w::Int) (str::[Int]) = + length (padTo 0 w str) >= length str+ && length (padTo 0 w str) >= w++quote str = '"':str++"\""++enumerate strs = zipWith (++) strs [show i | i<-[1..]]+numbered str = enumerate (repeat str)++wordsBy p s+ | findSpace == [] = []+ | otherwise = w : wordsBy p s''+ where+ (w, s'') = break p findSpace+ findSpace = dropWhile p s++printLabels ls = putStrLn $ concatMap (\(a,b) -> a++":\t"++show b++"\n") ls+
+ Utils/Table.hs view
@@ -0,0 +1,88 @@+module Utils.Table where+import Prelude hiding (map)+import Control.Monad+import qualified Data.List as L+import qualified Data.Map as Map+import Data.Function++type Table row col c = Map.Map (row,col) c+type Stripe a b = Map.Map a b++listTable lst = Map.fromList . concat $ [[((u,v),x) | x <- as | u <- [1..]] | as <- lst | v <- [1..]] ++fromList x = Map.fromList x+toList x = Map.toList x++elems = Map.elems+indexes = Map.keys++intersectionWith :: (Ord row, Ord col) => (a -> b -> c) + -> Table row col a+ -> Table row col b+ -> Table row col c+intersectionWith = Map.intersectionWith++-- map :: (c -> d) -> Table a b c -> Table a b d+map f = Map.map f +mapM f = seqTable . map f ++-- mapI :: ((a,b) -> c -> d) -> Table a b c -> Table a b d+mapI f = Map.mapWithKey f+mapIM f = seqTable . mapI f+mapIM_ f = seqTable_ . mapI f++row :: (Ord a, Ord b) => a -> Table a b c -> Stripe b c+row rowName table = rows table Map.! rowName ++column :: (Ord a, Ord b) => b -> Table a b c -> Stripe a c+column colName table = columns table Map.! colName ++rows :: (Ord a, Ord b) => Table a b c -> Map.Map a (Stripe b c)+rows = getStripe (snd,fst)++columns :: (Ord a, Ord b) => Table a b c -> Map.Map b (Stripe a c)+columns = getStripe (fst,snd)++rowNames :: (Ord a, Ord b) => Table a b c -> [a]+rowNames = L.nub . L.map fst . Map.keys++colNames :: (Ord a, Ord b) => Table a b c -> [b]+colNames = L.nub . L.map snd . Map.keys -- This is inefficient (nub)+++-- fromRows expects to form a full matrix+fromRows :: (Ord a, Ord b) => Map.Map a (Stripe b c) -> Table a b c+fromRows rs = Map.fromList+ [((r,c),rs Map.! r Map.! c) + | r <- Map.keys rs + , c <- L.nub $ L.concatMap Map.keys $ Map.elems rs]++fromCols :: (Ord a, Ord b) => Map.Map a (Stripe b c) -> Table b a c+fromCols rs = Map.fromList+ [((c,r),rs Map.! r Map.! c) + | r <- Map.keys rs + , c <- L.nub $ L.concatMap Map.keys $ Map.elems rs]+--++namedGroupsBy op = L.map (\x -> (op (head x),x))+ . L.groupBy (\x y -> op x == op y)+ . L.sortBy (compare `on` op)++getStripe (want,lose) table + = Map.fromList+ . L.map makeStripe+ . L.map dropR+ . namedGroupsBy (\(k,_) -> lose k)+ . Map.toList $ table+ where + dropR (n,s) = (n,L.map (\(k,v) -> (want k,v)) s)+ makeStripe (s,g) = (s,(Map.fromList g))+++-- seqTable :: (Ord a,Ord b,Monad m) => Table a b (m c) -> m (Table a b c)+seqTable table = Control.Monad.mapM (\(k,v) -> v >>= \r-> return (k,r)) (Map.toList table)+ >>= return . Map.fromList++seqTable_ table = Control.Monad.mapM_+ (\(k,v) -> v) (Map.toList table)+
+ Utils/Vector.hs view
@@ -0,0 +1,151 @@+module Utils.Vector where+import Data.List+import qualified Utils.List as UL+import Numeric+-- This is a quickie module for simple vector arithmetic on lists+-- It is provided so that optimization framework is not dependent+-- on GSL++type Vector = [Double]++-- The basic arithmetic+a <> b = sum $ zipWith (*) a b+s *| v = map (*s) v+a <-> b = zipWith (-) a b +a <+> b = zipWith (+) a b +a <*> b = zipWith (*) a b +v /| s = map (/s) v++infixl 7 <>, <*>, *|, /|+infixl 6 <+>, <->+++average vs = map ((/genericLength vs).sum) $ transpose vs+stdDev :: [Vector] -> Double+stdDev vs = UL.average $ map (\x -> norm (x <-> avg)) vs+ where avg = average vs++-- Projecting a vector to another+proj u v = ((v <> u) / (u <> u)) *| u++-- Norms and normalization+norm x = sqrt(x<>x)+normalize x = x /| norm x++normalizeN20 = normalize . map n . normalize+ where+ n x | x >0.2 = 0.2+ | otherwise = x ++-- Gram-Schmidt orthogonalization procedure+--gramSchmidt :: [[Double]] -> [[Double]] -> [[Double]]+gramSchmidt x = gs [] x+gs us [] = reverse us +gs [] (v:vs) = gs [v] vs+gs us (v:vs) = gs (vr:us) vs+ where vr = foldl1 (<->) $ v:[proj ui v | ui <- us]++-- Making an "identity matrix"+makeBase n = [[if y==x then 1 else 0 | x <- [1..n]] | y <- [1..n]]++-- a version of <= that works in different kind of wrong way than <=+safeLessThan a b = ( abs (a - b) <= 1.0e-9 ) || ( a < b) +++-- Utilities for printing vectors++showVector ::Int -> Vector -> String+showVector p l = "["++concat (intersperse "," (map (\n -> showFFloat (Just p) n "") l))++"]"++showMatrix1 :: Int -> [Vector] -> String+showMatrix1 p l = "["++concat (intersperse "," (map (showVector p) l))++"]"++showMatrix :: Int -> [Vector] -> String+showMatrix p l = "\t["++concat (intersperse ",\n\t" (map (showVector p) l))++"]"+++scaleToMax l v | n > l = l *| (v /| n)+ | otherwise = v+ where+ n = norm v+++-- Stuff dealing with ODE:s++rungeKutta4 h f (x0,t0) = (x0 <+> ((1/6) *| k1) + <+> ((1/3) *| k2)+ <+> ((1/3) *| k3)+ <+> ((1/6) *| k4)+ ,t0+h)+ where + k1 = h *| f x0 t0+ k2 = h *| f (x0<+> (0.5 *| k1)) (t0+h/2)+ k3 = h *| f (x0<+> (0.5 *| k2)) (t0+h/2)+ k4 = h *| f (x0<+> k3) (t0+h)++{-+class Vectorizeable a where+ toVector :: a -> Vector+ fromVector :: a -> Vector -> a++instance Vectorizeable Vector where+ toVector = id+ fromVector a = id+-}+{-+-- This is no good I think:+class ODE a b | a -> b where+ add :: a -> Double -> b -> a++addp :: (ODE a d) => a -> (Double,d) -> a+addp x (h,d) = add x h d++instance ODE Vector Vector where+ add a h d = a <+> (h*|d)+instance (ODE a b) => ODE [a] [b] where+ add a h d = zipWith (\ai di -> add ai h di) a d++eulerM h f (x0,t0) = do+ fx <- f x0 t0+ return (x0 <+> (h *| (f x0 t0)),t0+h)++eulerMODE :: (Monad m,ODE a d) => Double -> (a -> Double -> m d) + -> (a,Double) -> m (a,Double)+eulerMODE h f (x0,t0) = do+ fx <- f x0 t0+ return (x0 `addp` (h,fx),t0+h)++rungeKutta4ODE :: (ODE point delta) => Double -> (point -> Double -> delta) + -> (point,Double) -> (point,Double)+rungeKutta4ODE h f (x0,t0) = (x0 `addp` ((h/6),k1) + `addp` ((h/3),k2)+ `addp` ((h/3),k3)+ `addp` ((h/6),k4)+ ,t0+h)+ where + k1 = f x0 t0+ k2 = f (add x0 0.5 k1) (t0+h/2)+ k3 = f (add x0 0.5 k2) (t0+h/2)+ k4 = f (add x0 1 k3) (t0+h)++rungeKutta4MODE :: (Monad m, ODE point delta) => Double -> (point -> Double -> m delta) + -> (point,Double) -> m (point,Double)+rungeKutta4MODE h f (x0,t0) = do + k1 <- f x0 t0+ k2 <- f (add x0 0.5 k1) (t0+h/2)+ k3 <- f (add x0 0.5 k2) (t0+h/2)+ k4 <- f (add x0 1 k3) (t0+h)+ return (x0 `addp` ((h/6),k1) + `addp` ((h/3),k2)+ `addp` ((h/3),k3)+ `addp` ((h/6),k4)+ ,t0+h)+-}+-- Other numerical utilities. TODO: Split into own file++-- Snap a point to a grid+snap :: Double -> Double -> Double+snap s x = (fromIntegral $ round (x/is))*is where is = s++width (a,b) = abs (a-b)+