diff --git a/JYU-Utils.cabal b/JYU-Utils.cabal
new file mode 100644
--- /dev/null
+++ b/JYU-Utils.cabal
@@ -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
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/Utils/BinaryInstances.hs b/Utils/BinaryInstances.hs
new file mode 100644
--- /dev/null
+++ b/Utils/BinaryInstances.hs
@@ -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
diff --git a/Utils/Darcs.hs b/Utils/Darcs.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Darcs.hs
@@ -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"
diff --git a/Utils/DynMap.hs b/Utils/DynMap.hs
new file mode 100644
--- /dev/null
+++ b/Utils/DynMap.hs
@@ -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
+
+
diff --git a/Utils/Exception.hs b/Utils/Exception.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Exception.hs
@@ -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)) 
diff --git a/Utils/File.hs b/Utils/File.hs
new file mode 100644
--- /dev/null
+++ b/Utils/File.hs
@@ -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  "")
+
+
diff --git a/Utils/Function.hs b/Utils/Function.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Function.hs
@@ -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
diff --git a/Utils/List.hs b/Utils/List.hs
new file mode 100644
--- /dev/null
+++ b/Utils/List.hs
@@ -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
diff --git a/Utils/Monad.hs b/Utils/Monad.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Monad.hs
@@ -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
diff --git a/Utils/MonadRandom.hs b/Utils/MonadRandom.hs
new file mode 100644
--- /dev/null
+++ b/Utils/MonadRandom.hs
@@ -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
diff --git a/Utils/Numeric.hs b/Utils/Numeric.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Numeric.hs
@@ -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
+-}            
diff --git a/Utils/Parallel.hs b/Utils/Parallel.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Parallel.hs
@@ -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
+
+
diff --git a/Utils/Point.hs b/Utils/Point.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Point.hs
@@ -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)
diff --git a/Utils/Rectangle.hs b/Utils/Rectangle.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Rectangle.hs
@@ -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)
+
diff --git a/Utils/Sampling.hs b/Utils/Sampling.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Sampling.hs
@@ -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
diff --git a/Utils/SemanticEditors.hs b/Utils/SemanticEditors.hs
new file mode 100644
--- /dev/null
+++ b/Utils/SemanticEditors.hs
@@ -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
+
+
diff --git a/Utils/Shuffle.hs b/Utils/Shuffle.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Shuffle.hs
@@ -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)
diff --git a/Utils/ShuffleMWC.hs b/Utils/ShuffleMWC.hs
new file mode 100644
--- /dev/null
+++ b/Utils/ShuffleMWC.hs
@@ -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)
diff --git a/Utils/Stream.hs b/Utils/Stream.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Stream.hs
@@ -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
+			 
diff --git a/Utils/String.hs b/Utils/String.hs
new file mode 100644
--- /dev/null
+++ b/Utils/String.hs
@@ -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
+ 
diff --git a/Utils/Table.hs b/Utils/Table.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Table.hs
@@ -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)
+
diff --git a/Utils/Vector.hs b/Utils/Vector.hs
new file mode 100644
--- /dev/null
+++ b/Utils/Vector.hs
@@ -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)
+