diff --git a/forsyde-shallow.cabal b/forsyde-shallow.cabal
--- a/forsyde-shallow.cabal
+++ b/forsyde-shallow.cabal
@@ -1,5 +1,5 @@
 name:           forsyde-shallow
-version:        3.3.3.0
+version:        3.4.0.0
 cabal-version:  >= 1.8
 build-type:     Simple
 license:        BSD3
@@ -21,6 +21,7 @@
 tested-with:    GHC==7.10.3
               , GHC==8.0.2
               , GHC==8.2.1
+              , GHC==8.4.4
 
 -- In order to include all this files with sdist
 extra-source-files: LICENSE,
@@ -43,7 +44,6 @@
                  , ForSyDe.Shallow.Core.Signal
                  , ForSyDe.Shallow.Core.Vector
                  , ForSyDe.Shallow.Core.AbsentExt
-                 , ForSyDe.Shallow.Core.BitVector
                  , ForSyDe.Shallow.MoC
                  , ForSyDe.Shallow.MoC.CT
                  , ForSyDe.Shallow.MoC.Synchronous
@@ -66,14 +66,46 @@
                  , ForSyDe.Shallow.Utility.FilterLib
                  , ForSyDe.Shallow.Utility.Gaussian
                  , ForSyDe.Shallow.Utility.Memory
+                 , ForSyDe.Shallow.Utility.Matrix
+                 , ForSyDe.Shallow.Utility.BitVector
+  other-modules:
+      Paths_forsyde_shallow
   ghc-options:	-Wall -fno-warn-name-shadowing
 
 
 Test-Suite unit
   type:           exitcode-stdio-1.0
   hs-source-dirs: test
-  main-is:        Spec.hs
-  build-depends:  base>=4 && <6
+  main-is:        Unit.hs
+  build-depends:  base>=4.6 && <6
                 , forsyde-shallow
                 , hspec >= 2.2.4
   ghc-options:  -threaded -rtsopts -with-rtsopts=-N
+  other-modules:
+      Paths_forsyde_shallow
+
+test-suite doctests
+  type:             exitcode-stdio-1.0
+  hs-source-dirs:   test
+  main-is:          TestDocs.hs
+  ghc-options:      -threaded -isrc
+  build-depends:    base>=4.6 && <6
+                  , doctest >= 0.8
+                  , forsyde-shallow
+  other-modules:
+      Paths_forsyde_shallow
+
+test-suite properties
+  type:             exitcode-stdio-1.0
+  hs-source-dirs:   test
+  main-is:          Props.hs
+  ghc-options:      -threaded 
+  build-depends:    base>=4.6 && <6
+                  , QuickCheck
+                  , forsyde-shallow  
+                  , old-time
+                  , directory >= 1.2.2             
+                  , process >= 1.2.3   
+                  , random
+  other-modules:
+      Paths_forsyde_shallow
diff --git a/src/ForSyDe/Shallow/Core/AbsentExt.hs b/src/ForSyDe/Shallow/Core/AbsentExt.hs
--- a/src/ForSyDe/Shallow/Core/AbsentExt.hs
+++ b/src/ForSyDe/Shallow/Core/AbsentExt.hs
@@ -71,9 +71,8 @@
 
 psi = abstExtFunc
 
-
-
-
+instance Functor AbstExt where
+  fmap = abstExtFunc
 
 
 
diff --git a/src/ForSyDe/Shallow/Core/BitVector.hs b/src/ForSyDe/Shallow/Core/BitVector.hs
deleted file mode 100644
--- a/src/ForSyDe/Shallow/Core/BitVector.hs
+++ /dev/null
@@ -1,122 +0,0 @@
------------------------------------------------------------------------------
--- |
--- Module  :  ForSyDe.Shallow.Core.BitVector
--- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
--- License     :  BSD-style (see the file LICENSE)
--- 
--- Maintainer  :  forsyde-dev@ict.kth.se
--- Stability   :  experimental
--- Portability :  portable
---
--- It defines the bit vector operations from\/to integer.
------------------------------------------------------------------------------
-module ForSyDe.Shallow.Core.BitVector(
-    -- *Polynomial data type
-    BitVector, Parity(..),
-    -- *Bit-vector and integer transformations
-    intToBitVector, bitVectorToInt,
-    -- *Add even\/odd parity bit
-    addEvenParityBit, addOddParityBit, addParityBit,
-    -- *Remove parity bit
-    removeParityBit,
-    -- *Check the even\/odd parity of the bit-vector
-    isEvenParity, isOddParity,
-    -- *Judge the form sanity of the bit-vector
-    isBitVector
-  )
-  where
-
-import ForSyDe.Shallow.Core.Vector
-
-type BitVector = Vector Integer
-
--- |To judge whether the input bit-vector is in a proper form.
-isBitVector :: (Num t, Eq t) =>
-      Vector t  -- ^Input bit-vector
-   -> Bool  -- ^Output boolean value
-isBitVector NullV = True
-isBitVector (x:>xs) = (x `elem` [0, 1]) && isBitVector xs
-
--- |To transform the input integer to a bit-vector with specified number of
--- bits.
-intToBitVector :: Int  -- ^Num of the bits
-       -> Integer  -- ^The input integer
-       -> BitVector -- ^The output bit-vector
-intToBitVector bits n | n >= 0 && n < 2^(bits-1) 
-       = intToBitVector' bits n
-intToBitVector bits n | n < 0 && abs n <= 2^(bits-1) 
-       = intToBitVector' bits (n + 2^bits)
-intToBitVector _    _ | otherwise = 
-      error "intToBitvector : Number out of range!" 
-
--- |Helper function of 'intToBitVector'.
-intToBitVector' :: (Num a, Ord a1, Num a1, Integral t) => 
-       t -> a1 -> Vector a
-intToBitVector' 0    _ = NullV
-intToBitVector' bits n = if n >= 2^(bits-1) then
-         1 :> intToBitVector' (bits-1) (n - 2^(bits-1))
-       else  
-         0 :> intToBitVector' (bits-1) n
-
--- |To transform the input bit-vecotr to an integer.
-bitVectorToInt :: BitVector -> Integer
-bitVectorToInt (1:>xv) | isBitVector xv 
-       = bitVectorToInt' xv (lengthV xv) - 2 ^ lengthV xv 
-bitVectorToInt (0:>xv) | isBitVector xv 
-       = bitVectorToInt' xv (lengthV xv)
-bitVectorToInt _ = error "bitVectorToInt: Vector is not a BitVector!"
-
-
--- |Helper function of 'bitVectorToInt'.
-bitVectorToInt' :: (Integral a, Num t) => Vector t -> a -> t
-bitVectorToInt' NullV   _   = 0
-bitVectorToInt' (x:>xv) bit = x * 2^(bit-1) + bitVectorToInt' xv (bit-1)
-
-data Parity = Even | Odd deriving (Show, Eq)
-
--- |To add even parity bit on the bit-vector in the tail.
-addEvenParityBit :: (Num a, Eq a) => Vector a -> Vector a
-addEvenParityBit = addParityBit Even
--- |To add odd parity bit on the bit-vector in the tail.
-addOddParityBit :: (Num a, Eq a) => Vector a -> Vector a
-addOddParityBit  = addParityBit Odd
-
-addParityBit :: (Num a, Eq a) => Parity -> Vector a -> Vector a
-addParityBit p v 
-  | isBitVector v = case p of
-    Even -> resZero even
-    Odd -> resZero (not even)
-  | otherwise =  error "addParity: Vector is not a BitVector"  
- where even = evenNumber v 
-       resZero b = v <+> unitV (if b then 0 else 1)
-
-
--- |To remove the parity bit in the tail.
-removeParityBit :: (Num t, Eq t) => Vector t -> Vector t
-removeParityBit v 
- | isBitVector v = takeV (lengthV v - 1) v
- | otherwise = error "removeParityBit: Vector is not a BitVector "
-
--- |To check the even parity of the bit-vector.
-isEvenParity :: (Num t, Eq t) => Vector t -> Bool
-isEvenParity = isParityCorrect Even
-
--- |To check the odd parity of the bit-vector.
-isOddParity :: (Num t, Eq t) => Vector t -> Bool
-isOddParity = isParityCorrect Odd
-
-isParityCorrect :: (Num t, Eq t) => Parity -> Vector t -> Bool 
-isParityCorrect Even xv = evenNumber xv
-isParityCorrect Odd xv  = not $ evenNumber xv 
-
-evenNumber :: (Num t, Eq t) => Vector t -> Bool
-evenNumber NullV   = True
-evenNumber (0:>xv) = xor False (evenNumber xv)
-evenNumber (1:>xv) = xor True (evenNumber xv)
-evenNumber (_:>_) = error "evenNumber: Vector is not a BitVector "
-         
-xor :: Bool -> Bool -> Bool
-xor True  False = True
-xor False True  = True
-xor _     _     = False 
-
diff --git a/src/ForSyDe/Shallow/Core/Vector.hs b/src/ForSyDe/Shallow/Core/Vector.hs
--- a/src/ForSyDe/Shallow/Core/Vector.hs
+++ b/src/ForSyDe/Shallow/Core/Vector.hs
@@ -1,140 +1,452 @@
 -----------------------------------------------------------------------------
 -- |
--- Module  :  ForSyDe.Shallow.Core.Vector
--- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
+-- Module      :  ForSyDe.Shallow.Core.Vector
+-- Copyright   :  (c) ForSyDe Group, KTH 2007-2019
 -- License     :  BSD-style (see the file LICENSE)
 -- 
--- Maintainer  :  forsyde-dev@ict.kth.se
+-- Maintainer  :  forsyde-dev@kth.se
 -- Stability   :  experimental
 -- Portability :  portable
 --
--- This module defines the data type 'Vector' and the
--- corresponding functions. It is a development of the module
--- defined by Reekie.  Though the vector is modeled as a list, it
--- should be viewed as an array, i.e. a vector has a fixed
--- size. Unfortunately, it is not possible to have the size of the
--- vector as a parameter of the vector data type, due to restrictions
--- in Haskells type system. Still most operations are defined for
--- vectors with the same size.
+-- This module defines the data type 'Vector' and the corresponding
+-- functions. It is a development of the module defined by
+-- <https://ptolemy.berkeley.edu/~johnr/papers/pdf/thesis.pdf Reekie>.
+-- The 'Vector' data type is a shallow interpretation of arrays and is
+-- used for quick prototyping of array algorithms and skeletons,
+-- whereas in fact it is implemented as an infinite list itself. For a
+-- type-checked fixed-size data type for representing vectors, see
+-- <http://hackage.haskell.org/package/parameterized-data FSVec> or
+-- <http://hackage.haskell.org/package/repa REPA>.
+--
+-- __OBS:__ The lengths in the API documentation for function arguments
+-- are not type-safe, but rather suggestions for usage in designing
+-- vector algorithms or skeletons.
 -----------------------------------------------------------------------------
-module ForSyDe.Shallow.Core.Vector ( 
-  Vector (..), vector, fromVector, unitV, nullV, lengthV,
-  atV, replaceV, headV, tailV, lastV, initV, takeV, dropV, 
-  selectV, groupV, (<+>), (<:), mapV, foldlV, foldrV,
-  reduceV, pipeV, zipWithV, filterV, zipV, unzipV, 
-  -- scanlV, scanrV, meshlV, meshrV, 
-  concatV, reverseV, shiftlV, shiftrV, rotrV, rotlV, rotateV,
-  generateV, iterateV, copyV --, serialV, parallelV 
+module ForSyDe.Shallow.Core.Vector (
+  Vector (..), (<+>), (<:), 
+  -- * Queries
+  nullV, lengthV,
+  -- * Generators
+  vector, fromVector, unitV, 
+  iterateV, generateV, copyV,
+  -- * Functional skeletons
+  mapV, zipWithV, zipWith3V,
+  reduceV, pipeV, foldlV, foldrV, 
+  scanlV, scanrV, -- meshlV, meshrV,
+  -- * Selectors
+  atV, headV, tailV, lastV, initV, headsV, tailsV,
+  takeV, dropV, selectV, groupV, filterV, stencilV,
+  -- * Permutators
+  replaceV, zipV, unzipV,
+  concatV, reverseV, shiftlV, shiftrV, rotrV, rotlV, rotateV
   ) where
 
+-----------------------------------------------------------------------------
+-- CONSTRUCTORS AND INSTANCES
+-----------------------------------------------------------------------------
+
 infixr 5 :>
 infixl 5 <:
 infixr 5 <+>
 
--- | The data type 'Vector' is modeled similar to a list. It has two data type constructors. 'NullV' constructs the empty vector, while ':>' constructsa vector by adding an value to an existing vector. Using the inheritance mechanism of Haskell we have declared 'Vector' as an instance of the classes 'Read' and 'Show'.
+-- | The data type 'Vector' is modeled similar to a list. It has two data type constructors. 'NullV' constructs the empty vector, while ':>' constructsa vector by adding an value to an existing vector..
 --
--- | This means that the vector 1:>2:>3:>NullV is shown as <1,2,3>.
+-- 'Vector' is an instance of the classes 'Read' and 'Show'. This means that the vector
+--
+-- > 1:>2:>3:>NullV
+--
+-- is shown as
+--
+-- > <1,2,3>
 data Vector a = NullV
               | a :> (Vector a) deriving (Eq)
 
+instance (Show a) => Show (Vector a) where
+  showsPrec p NullV = showParen (p > 9) (showString "<>")
+  showsPrec p xs    = showParen (p > 9) (showChar '<' . showVector1 xs)
+    where
+      showVector1 NullV = showChar '>'            
+      showVector1 (y:>NullV) = shows y . showChar '>'
+      showVector1 (y:>ys) = shows y . showChar ',' 
+                            . showVector1 ys
+
+instance Read a => Read (Vector a) where
+   readsPrec _ s = readsVector s
+
+readsVector :: (Read a) => ReadS (Vector a)
+readsVector s = [((x:>NullV), rest) | ("<", r2) <- lex s,
+            (x, r3)   <- reads r2,
+            (">", rest) <- lex r3]
+       ++
+      [(NullV, r4)    | ("<", r5) <- lex s,
+            (">", r4) <- lex r5]
+       ++
+      [((x:>xs), r6)  | ("<", r7) <- lex s,
+            (x, r8)   <- reads r7,
+            (",", r9) <- lex r8,
+            (xs, r6) <- readsValues r9]
+
+readsValues :: (Read a) => ReadS (Vector a)
+readsValues s = [((x:>NullV), r1) | (x, r2)   <- reads s,
+              (">", r1) <- lex r2]
+      ++
+      [((x:>xs), r3)    | (x, r4)   <- reads s,
+              (",", r5) <- lex r4,
+              (xs, r3)  <- readsValues r5]
+
+-- | The operator '(<:)' appends an element at the end of a vector.
+(<:)  :: Vector a  -- ^ /length/ = @la@
+      -> a
+      -> Vector a  -- ^ /length/ = @la + 1@
+xs <: x = xs <+> unitV x     
+
+-- | The operator '<+>' concatenates two vectors.
+(<+>) :: Vector a  -- ^ /length/ = @la@
+      -> Vector a  -- ^ /length/ = @lb@
+      -> Vector a  -- ^ /length/ = @la + lb@
+NullV <+> ys   = ys
+(x:>xs) <+> ys = x :> (xs <+> ys) 
+
+-----------------------------------------------------------------------------
+-- GENERATORS
+-----------------------------------------------------------------------------
+
 -- | The function 'vector' converts a list into a vector.
-vector     :: [a] -> Vector a
+vector        :: [a] -> Vector a
+vector []     = NullV
+vector (x:xs) = x :> (vector xs)
 
 -- | The function 'fromVector' converts a vector into a list.
-fromVector :: Vector a -> [a]
+fromVector         :: Vector a -> [a]
+fromVector NullV   = []
+fromVector (x:>xs) = x : fromVector xs
 
 -- | The function 'unitV' creates a vector with one element. 
-unitV :: a -> Vector a
+unitV   :: a -> Vector a  -- ^ /length/ = @1@
+unitV x = x :> NullV
 
+-- | The function 'iterateV' generates a vector with a given number of
+-- elements starting from an initial element using a supplied function
+-- for the generation of elements.
+--
+-- >>> iterateV 5 (+1) 1
+-- <1,2,3,4,5>
+iterateV :: (Num a, Eq a)
+         => a        -- ^ number of elements = @n@
+         -> (b -> b) -- ^ generator function (@last_element -> next_element@)
+         -> b        -- ^ initial element
+         -> Vector b -- ^ generated vector; /length/ = @n@
+iterateV 0 _ _ = NullV
+iterateV n f a = a :> iterateV (n-1) f (f a)
+
+-- | The function 'generateV' behaves in the same way as 'iterateV',
+-- but starts with the application of the supplied function to the
+-- supplied value.
+--
+-- >>> generateV 5 (+1) 1
+-- <2,3,4,5,6>
+generateV :: (Num a, Eq a)
+         => a        -- ^ number of elements = @n@
+         -> (b -> b) -- ^ generator function (@last_element -> next_element@)
+         -> b        -- ^ initial element
+         -> Vector b -- ^ generated vector; /length/ = @n@
+generateV 0 _ _ = NullV
+generateV n f a = x :> generateV (n-1) f x 
+        where x = f a
+
+-- | The function 'copyV' generates a vector with a given number of
+-- copies of the same element.
+--
+-- >>> copyV 7 5 
+-- <5,5,5,5,5,5,5>
+copyV     :: (Num a, Eq a)
+          => a        -- ^ number of elements = @n@
+          -> b        -- ^ element to be copied
+          -> Vector b -- ^ /length/ = @n@
+copyV k x = iterateV k id x 
+
+-----------------------------------------------------------------------------
+-- QUERIES
+-----------------------------------------------------------------------------
+
 -- | The function 'nullV' returns 'True' if a vector is empty. 
-nullV :: Vector a -> Bool
+nullV       :: Vector a -> Bool
+nullV NullV = True
+nullV _     = False
 
 -- | The function 'lengthV' returns the number of elements in a value. 
-lengthV :: Vector a -> Int
+lengthV         :: Vector a -> Int
+lengthV NullV   = 0
+lengthV (_:>xs) = 1 + lengthV xs
 
--- | The function 'atV' returns the n-th element in a vector, starting from zero.
-atV  :: (Num a, Eq a) => Vector b -> a -> b
+-----------------------------------------------------------------------------
+-- HIGHER ORDER SKELETONS
+-----------------------------------------------------------------------------
 
--- |  The function 'replaceV' replaces an element in a vector.
-replaceV :: Vector a -> Int -> a -> Vector a
+-- | The higher-order function 'mapV' applies a function on all elements of a vector.
+mapV :: (a -> b)
+     -> Vector a  -- ^ /length/ = @la@
+     -> Vector b  -- ^ /length/ = @la@
+mapV f (x:>xs) = f x :> mapV f xs
+mapV _ NullV   = NullV
 
+-- | The higher-order function 'zipWithV' applies a function pairwise on two vectors.
+zipWithV :: (a -> b -> c)
+         -> Vector a  -- ^ /length/ = @la@
+         -> Vector b  -- ^ /length/ = @lb@
+         -> Vector c  -- ^ /length/ = @minimum [la,lb]@
+zipWithV f (x:>xs) (y:>ys) = f x y :> (zipWithV f xs ys)
+zipWithV _ _ _ = NullV
+
+-- | The higher-order function 'zipWithV3' applies a function 3-tuple-wise on three vectors.
+zipWith3V :: (a -> b -> c -> d)
+          -> Vector a  -- ^ /length/ = @la@
+          -> Vector b  -- ^ /length/ = @lb@
+          -> Vector c  -- ^ /length/ = @lc@
+          -> Vector d  -- ^ /length/ = @minimum [la,lb,lc]@
+zipWith3V f (x:>xs) (y:>ys) (z:>zs) = f x y z :> (zipWith3V f xs ys zs)
+zipWith3V _ _ _ _ = NullV
+
+-- | The higher-order functions 'foldlV' folds a function from the
+-- right to the left over a vector using an initial value.
+--
+-- >>> foldlV (-) 8 $ vector [4,2,1]   -- is the same as (((8 - 4) - 2) - 1) 
+-- 1
+foldlV :: (a -> b -> a) -> a -> Vector b -> a 
+foldlV _ a NullV   = a
+foldlV f a (x:>xs) = foldlV f (f a x) xs
+
+-- | The higher-order functions 'foldrV' folds a function from the
+-- left to the right over a vector using an initial value.
+--
+-- >>> foldrV (-) 8 $ vector [4,2,1]   -- is the same as (4 - (2 - (1 - 8)))
+-- -5
+foldrV :: (b -> a -> a) -> a -> Vector b -> a
+foldrV _ a NullV   = a 
+foldrV f a (x:>xs) = f x (foldrV f a xs)
+
+-- | Reduces a vector of elements to a single element based on a
+-- binary function.
+--
+-- >>> reduceV (+) $ vector [1,2,3,4,5]
+-- 15
+reduceV :: (a -> a -> a) -> Vector a -> a
+reduceV _ NullV      = error "Cannot reduce a null vector"
+reduceV _ (x:>NullV) = x
+reduceV f (x:>xs)    = foldlV f x xs
+
+-- | Pipes an element through a vector of functions.
+--
+-- >>> vector [(*2), (+1), (/3)] `pipeV` 3      -- is the same as ((*2) . (+1) . (/3)) 3
+-- 4.0
+pipeV :: Vector (a -> a) -> a -> a
+pipeV vf = foldrV (.) id vf
+
+-----------------------------------------------------------------------------
+-- SELECTORS
+-----------------------------------------------------------------------------
+
+-- | The function 'atV' returns the n-th element in a vector, starting
+-- from zero.
+--
+-- >>> vector [1,2,3,4,5] `atV` 3
+-- 4
+atV  :: (Integral a) => Vector b -> a -> b
+NullV   `atV` _ = error "atV: Vector has not enough elements"
+(x:>_)  `atV` 0 = x
+(_:>xs) `atV` n = xs `atV` (n-1)
+
 -- | The functions 'headV' returns the first element of a vector.
 headV :: Vector a -> a
+headV NullV   = error "headV: Vector is empty"
+headV (v:>_) = v
 
 -- | The function 'lastV' returns the last element of a vector.
 lastV :: Vector a -> a
+lastV NullV  = error "lastV: Vector is empty"
+lastV (v:>NullV) = v
+lastV (_:>vs)    = lastV vs
 
 -- | The functions 'tailV' returns all, but the first element of a vector.
-tailV :: Vector a -> Vector a 
+tailV :: Vector a  -- ^ /length/ = @la@
+      -> Vector a  -- ^ /length/ = @la-1@
+tailV NullV   = error "tailV: Vector is empty"
+tailV (_:>vs) = vs
 
 -- | The function 'initV' returns all but the last elements of a vector.
-initV :: Vector a -> Vector a 
-
--- | The function 'takeV' returns the first n elements of a vector.
-takeV :: (Num a, Ord a) => a -> Vector b -> Vector b
-
--- | The function 'dropV' drops the first n elements of a vector.
-dropV :: (Num a, Ord a) => a -> Vector b -> Vector b
+initV :: Vector a  -- ^ /length/ = @la@
+      -> Vector a  -- ^ /length/ = @la-1@
+initV NullV  = error "initV: Vector is empty"
+initV (_:>NullV) = NullV
+initV (v:>vs)    = v :> initV vs
 
--- | The function 'selectV' selects elements in the vector. The first argument gives the initial element, starting from zero, the second argument gives the stepsize between elements and the last argument gives the number of elements. 
-selectV :: Int -> Int -> Int -> Vector a -> Vector a
+-- | The function 'takeV' returns the first @n@ elements of a vector.
+-- 
+-- >>> takeV 2 $ vector [1,2,3,4,5]
+-- <1,2>
+takeV :: (Num a, Ord a)
+      => a        -- ^ @= n@
+      -> Vector b -- ^ /length/ = @la@
+      -> Vector b -- ^ /length/ = @minimum [n,la]@
+takeV 0 _       = NullV
+takeV _ NullV       = NullV
+takeV n (v:>vs) | n <= 0    = NullV
+                | otherwise = v :> takeV (n-1) vs
 
--- | The function 'groupV' groups a vector into a vector of vectors of size n.
-groupV :: Int -> Vector a -> Vector (Vector a)
+-- | The function 'dropV' drops the first @n@ elements of a vector.
+--
+-- >>> dropV 2 $ vector [1,2,3,4,5]
+-- <3,4,5>
+dropV :: (Num a, Ord a)
+      => a        -- ^ @= n@
+      -> Vector b -- ^ /length/ = @la@
+      -> Vector b -- ^ /length/ = @maximum [0,la-n]@
+dropV 0 vs      = vs
+dropV _ NullV       = NullV
+dropV n (v:>vs) | n <= 0    = v :> vs
+                | otherwise = dropV (n-1) vs
 
--- | The operator '(<:)' adds an element at the end of a vector.
-(<:)  :: Vector a -> a -> Vector a
+-- | The function 'selectV' selects elements in the vector based on a
+-- regular stride.
+selectV :: Int      -- ^ the initial element, starting from zero
+        -> Int      -- ^ stepsize between elements
+        -> Int      -- ^ number of elements @= n@
+        -> Vector a -- ^ /length/ = @la@ 
+        -> Vector a -- ^ /length/ @= n@
+selectV f s n vs
+  | n <= 0                 = NullV
+  | (f+s*n-1) > lengthV vs = error "selectV: Vector has not enough elements"
+  | otherwise              = atV vs f :> selectV (f+s) s (n-1) vs
 
--- | The operator '(<+>)' concatinates two vectors.
-(<+>) :: Vector a -> Vector a -> Vector a
+-- | The function 'groupV' groups a vector into a vector of vectors of
+-- size n.
+--
+-- >>> groupV 3 $ vector [1,2,3,4,5,6,7,8]
+-- <<1,2,3>,<4,5,6>>
+groupV :: Int               -- ^ @= n@
+       -> Vector a          -- ^ /length/ = @la@ 
+       -> Vector (Vector a) -- ^ /length/ = @la `div` n@ 
+groupV n v 
+  | lengthV v < n = NullV
+  | otherwise     = selectV 0 1 n v 
+                    :> groupV n (selectV n 1 (lengthV v-n) v)
 
 
--- | The higher-order function 'mapV' applies a function on all elements of a vector.
-mapV :: (a -> b) -> Vector a -> Vector b    
+-- | The higher-function 'filterV' takes a predicate function and a
+-- vector and creates a new vector with the elements for which the
+-- predicate is true.
+--
+-- >>> filterV odd $ vector [1,2,3,4,5,6,7,8]
+-- <1,3,5,7>
+--
+-- (*) however, the length is __unknown__, because it is dependent on
+-- the data contained inside the vector. Try avoiding 'filterV' in
+-- designs where the size of the data is crucial.
+filterV :: (a -> Bool) -- ^ predicate function
+        -> Vector a    -- ^ /length/ = @la@
+        -> Vector a    -- ^ /length/ @<= la@ (*)
+filterV _ NullV   = NullV
+filterV p (v:>vs) = if (p v)
+                    then v :> filterV p vs
+                    else filterV p vs
 
--- | The higher-order function 'zipWithV' applies a function pairwise on to vectors.
-zipWithV :: (a -> b -> c) -> Vector a -> Vector b -> Vector c
+-- | Returns a vector containing all the possible prefixes of an input
+-- vector.
+--
+-- >>> let v = vector [1,2,3,4,5,6]
+-- >>> headsV v
+-- <<1>,<1,2>,<1,2,3>,<1,2,3,4>,<1,2,3,4,5>,<1,2,3,4,5,6>,<1,2,3,4,5,6>>
+headsV :: Vector a          -- ^ /length/ = @la@
+       -> Vector (Vector a) -- ^ /length/ = @la + 1@
+headsV NullV  = error "heads: null vector"
+headsV v      = foldrV sel (unitV NullV) $ mapV (unitV . unitV) v
+  where sel x y = x <+> mapV (lastV  x <+>) y
 
--- | The higher-order functions 'foldlV' folds a function from the right to the left  over a vector using an initial value.
-foldlV :: (a -> b -> a) -> a -> Vector b -> a 
+-- | Returns a vector containing all the possible suffixes of an input
+-- vector.
+--
+-- >>> let v = vector [1,2,3,4,5,6]
+-- >>> tailsV v
+-- <<1,2,3,4,5,6>,<2,3,4,5,6>,<3,4,5,6>,<4,5,6>,<5,6>,<6>,<>>
+tailsV :: Vector a          -- ^ /length/ = @la@
+       -> Vector (Vector a) -- ^ /length/ = @la + 1@
+tailsV NullV = NullV
+tailsV v    = foldrV sel (unitV NullV) $ mapV (unitV . unitV) v
+  where sel x y = mapV (<+> headV y) x <+> y
 
--- | The higher-order functions 'foldrV' folds a function from the left to the right over a vector using an initial value.
-foldrV :: (b -> a -> a) -> a -> Vector b -> a
+-- | Returns a stencil of @n@ neighboring elements for each possible
+-- element in a vector.
+--
+-- >>> stencilV 3 $ vector [1..5]
+-- <<1,2,3>,<2,3,4>,<3,4,5>>
+stencilV :: Int               -- ^ stencil size @= n@
+         -> Vector a          -- ^ /length/ = @la@ 
+         -> Vector (Vector a) -- ^ /length/ = @la - n + 1@ 
+stencilV n v = mapV (takeV n) $ dropFromEnd n $ tailsV v
+  where dropFromEnd n = takeV (lengthV v - n + 1)
 
--- | Reduces a vector of elements to a single element based on a binary function. 
-reduceV :: (a -> a -> a) -> Vector a -> a
+-----------------------------------------------------------------------------
+-- PERMUTATORS
+-----------------------------------------------------------------------------
 
--- | Pipes an element through a vector of functions. For example the code
---
--- >>> pipeV [(*2), (+1), (/3)] 3
--- > 2
---
--- is equivalent to
+-- |  The function 'replaceV' replaces an element in a vector.
 --
--- >>> ((*2) . (+1) . (/3)) 3
--- > 2
-pipeV :: Vector (a -> a) -> a -> a
-
--- | The higher-function 'filterV' takes a predicate function and a vector and creates a new vector with the elements for which the predicate is true. 
-filterV :: (a -> Bool) -> Vector a -> Vector a
+-- >>> replaceV (vector [1..5]) 2 100
+-- <1,2,100,4,5>
+replaceV :: Vector a -- ^ input vector; /length/ = @la@
+         -> Int      -- ^ position of the element to be replaced
+         -> a        -- ^ new element
+         -> Vector a -- ^ altered vector; /length/ = @la@
+replaceV vs n x 
+    | n <= lengthV vs && n >= 0 = takeV n vs <+> unitV x 
+                                  <+> dropV (n+1) vs
+    | otherwise                 = vs
 
 -- | The function 'zipV' zips two vectors into a vector of tuples.
-zipV   :: Vector a -> Vector b -> Vector (a, b)
+zipV   :: Vector a      -- ^ /length/ = @la@ 
+       -> Vector b      -- ^ /length/ = @lb@ 
+       -> Vector (a, b) -- ^ /length/ = @minimum [la,lb]@ 
+zipV (x:>xs) (y:>ys) = (x, y) :> zipV xs ys
+zipV _   _   = NullV
 
 -- | The function 'unzipV' unzips a vector of tuples into two vectors.
-unzipV :: Vector (a, b) -> (Vector a, Vector b)
+unzipV :: Vector (a, b)        -- ^ /length/ = @la@
+       -> (Vector a, Vector b) -- ^ /length/ = @la@
+unzipV NullV           = (NullV, NullV)
+unzipV ((x, y) :> xys) = (x:>xs, y:>ys) 
+  where (xs, ys) = unzipV xys
 
--- | The function 'shiftlV' shifts a value from the left into a vector. 
-shiftlV :: Vector a -> a-> Vector a 
+-- | The function 'shiftlV' shifts a value from the left into a vector.
+--
+-- >>> vector [1..5] `shiftlV` 100
+-- <100,1,2,3,4>
+shiftlV :: Vector a -> a -> Vector a 
+shiftlV vs v = v :> initV vs
 
 -- | The function 'shiftrV' shifts a value from the right into a vector. 
+--
+-- >>> vector [1..5] `shiftrV` 100
+-- <2,3,4,5,100>
 shiftrV :: Vector a -> a -> Vector a
+shiftrV vs v = tailV vs <: v
 
--- | The function 'rotlV' rotates a vector to the left. Note that this fuctions does not change the size of a vector.
+-- | The function 'rotlV' rotates a vector to the left. Note that this
+-- fuctions does not change the size of a vector.
+--
+-- >>> rotlV $ vector [1..5]
+-- <5,1,2,3,4>
 rotlV   :: Vector a -> Vector a
+rotrV NullV = NullV
+rotrV vs    = tailV vs <: headV vs
 
--- | The function 'rotrV' rotates a vector to the right. Note that this fuction does not change the size of a vector.
+-- | The function 'rotrV' rotates a vector to the right. Note that
+-- this fuction does not change the size of a vector.
+--
+-- >>> rotrV $ vector [1..5]
+-- <2,3,4,5,1>
 rotrV   :: Vector a -> Vector a
+rotlV NullV = NullV
+rotlV vs    = lastV vs :> initV vs
 
 -- | The function 'rotateV' rotates a vector based on an index offset.
 --
@@ -146,33 +458,44 @@
 -- * @(< 0)@ : rotates the vector right with the corresponding number
 -- of positions.
 rotateV :: Int -> Vector a -> Vector a
+rotateV n
+  | n > 0     = pipeV (copyV (abs n) rotlV)
+  | n < 0     = pipeV (copyV (abs n) rotrV)
+  | otherwise = id
 
 -- | The function 'concatV' transforms a vector of vectors to a single vector. 
 concatV   :: Vector (Vector a) -> Vector a
+concatV = foldrV (<+>) NullV
 
 -- | The function 'reverseV' reverses the order of elements in a vector. 
 reverseV  :: Vector a -> Vector a
-
--- | The function 'iterateV' generates a vector with a given number of elements starting from an initial element using a supplied function for the generation of elements. 
---
--- > Vector> iterateV 5 (+1) 1
---
--- > <1,2,3,4,5> :: Vector Integer
-iterateV :: (Num a, Eq a) => a -> (b -> b) -> b -> Vector b
+reverseV NullV   = NullV
+reverseV (v:>vs) = reverseV vs <: v
 
--- | The function 'generateV' behaves in the same way, but starts with the application of the supplied function to the supplied value. 
+-- | Performs the parallel prefix operation on a vector.
 --
--- > Vector> generateV 5 (+1) 1
--- 
--- > <2,3,4,5,6> :: Vector Integer
-generateV :: (Num a, Eq a) => a -> (b -> b) -> b -> Vector b
-
--- | The function 'copyV' generates a vector with a given number of copies of the same element. 
+-- >>> scanlV (+) 0 $ vector [1,1,1,1,1,1]
+-- <1,2,3,4,5,6>
+scanlV    :: (a -> b -> a)  -- ^ funtion to generate next element
+          -> a              -- ^ initial element
+          -> Vector b       -- ^ input vector; /length/ = @l@
+          -> Vector a       -- ^ output vector; /length/ = @l@ 
+scanlV _ _ NullV   = NullV
+scanlV f a (x:>xs) = q :> scanlV f q xs 
+       where q = f a x
+             
+-- | Performs the parallel suffix operation on a vector.
 --
--- > Vector> copyV 7 5 
--- 
--- > <5,5,5,5,5,5,5> :: Vector Integer
-copyV     :: (Num a, Eq a) => a -> b -> Vector b
+-- >>> scanrV (+) 0 $ vector [1,1,1,1,1,1]
+-- <6,5,4,3,2,1>
+scanrV    :: (b -> a -> a)   -- ^ funtion to generate next element
+          -> a               -- ^ initial element       
+          -> Vector b        -- ^ input vector; /length/ = @l@
+          -> Vector a        -- ^ output vector; /length/ = @l@ 
+scanrV _ _ NullV  = NullV
+scanrV f a (x:>NullV) = f x a :> NullV
+scanrV f a (x:>xs)    = f x y :> ys 
+          where ys@(y:>_) = scanrV f a xs
 
 {-
 -- | The function 'serialV' can be used to construct a serial network of processes.
@@ -217,170 +540,7 @@
 \index{meshrV@\haskell{meshrV}}
 -}
 
-instance (Show a) => Show (Vector a) where
-   showsPrec p NullV = showParen (p > 9) (
-             showString "<>")
-   showsPrec p xs    = showParen (p > 9) (
-             showChar '<' . showVector1 xs)
-           where
-          showVector1 NullV
-             = showChar '>'            
-          showVector1 (y:>NullV) 
-             = shows y . showChar '>'
-          showVector1 (y:>ys)    
-             = shows y . showChar ',' 
-           . showVector1 ys
 
-
-instance Read a => Read (Vector a) where
-   readsPrec _ s = readsVector s
-
-readsVector :: (Read a) => ReadS (Vector a)
-readsVector s = [((x:>NullV), rest) | ("<", r2) <- lex s,
-            (x, r3)   <- reads r2,
-            (">", rest) <- lex r3]
-       ++
-      [(NullV, r4)    | ("<", r5) <- lex s,
-            (">", r4) <- lex r5]
-       ++
-      [((x:>xs), r6)  | ("<", r7) <- lex s,
-            (x, r8)   <- reads r7,
-            (",", r9) <- lex r8,
-            (xs, r6) <- readsValues r9]
-
-readsValues :: (Read a) => ReadS (Vector a)
-readsValues s = [((x:>NullV), r1) | (x, r2)   <- reads s,
-              (">", r1) <- lex r2]
-      ++
-      [((x:>xs), r3)    | (x, r4)   <- reads s,
-              (",", r5) <- lex r4,
-              (xs, r3)  <- readsValues r5]
-
-vector []     = NullV
-vector (x:xs) = x :> (vector xs)
-
-fromVector NullV   = []
-fromVector (x:>xs) = x : fromVector xs
-
-unitV x = x :> NullV
-
-nullV NullV   = True
-nullV _     = False
-
-lengthV NullV   = 0
-lengthV (_:>xs) = 1 + lengthV xs
-
-replaceV vs n x 
-    | n <= lengthV vs && n >= 0 = takeV n vs <+> unitV x 
-          <+> dropV (n+1) vs
-    | otherwise         =  vs
-
-NullV   `atV` _ = error "atV: Vector has not enough elements"
-(x:>_)  `atV` 0 = x
-(_:>xs) `atV` n = xs `atV` (n-1)
-
-headV NullV   = error "headV: Vector is empty"
-headV (v:>_) = v
-
-tailV NullV   = error "tailV: Vector is empty"
-tailV (_:>vs) = vs
-
-lastV NullV  = error "lastV: Vector is empty"
-lastV (v:>NullV) = v
-lastV (_:>vs)    = lastV vs
-
-initV NullV  = error "initV: Vector is empty"
-initV (_:>NullV) = NullV
-initV (v:>vs)    = v :> initV vs
-
-takeV 0 _       = NullV
-takeV _ NullV       = NullV
-takeV n (v:>vs) | n <= 0    = NullV
-        | otherwise = v :> takeV (n-1) vs
-
-dropV 0 vs      = vs
-dropV _ NullV       = NullV
-dropV n (v:>vs) | n <= 0    = v :> vs
-        | otherwise = dropV (n-1) vs
-
-selectV f s n vs | n <= 0       
-         = NullV
-         | (f+s*n-1) > lengthV vs 
-        = error "selectV: Vector has not enough elements"
-         | otherwise    
-        = atV vs f :> selectV (f+s) s (n-1) vs
-
-groupV n v 
-  | lengthV v < n = NullV
-  | otherwise     = selectV 0 1 n v 
-      :> groupV n (selectV n 1 (lengthV v-n) v)
-
-NullV <+> ys   = ys
-(x:>xs) <+> ys = x :> (xs <+> ys) 
-
-xs <: x = xs <+> unitV x     
-
-mapV _ NullV   = NullV
-mapV f (x:>xs) = f x :> mapV f xs
-
-zipWithV f (x:>xs) (y:>ys) = f x y :> (zipWithV f xs ys)
-zipWithV _ _   _   = NullV
-
-foldlV _ a NullV   = a
-foldlV f a (x:>xs) = foldlV f (f a x) xs
-
-foldrV _ a NullV   = a 
-foldrV f a (x:>xs) = f x (foldrV f a xs)
-
-reduceV _ NullV      = error "Cannot reduce a null vector"
-reduceV _ (x:>NullV) = x
-reduceV f (x:>xs)    = foldrV f x xs
-
-pipeV = reduceV (.)
-
-filterV _ NullV   = NullV
-filterV p (v:>vs) = if (p v) then
-         v :> filterV p vs
-      else 
-         filterV p vs
-
-zipV (x:>xs) (y:>ys) = (x, y) :> zipV xs ys
-zipV _   _   = NullV
-
-unzipV NullV       = (NullV, NullV)
-unzipV ((x, y) :> xys) = (x:>xs, y:>ys) 
-       where (xs, ys) = unzipV xys
-
-shiftlV vs v = v :> initV vs
-
-shiftrV vs v = tailV vs <: v
-
-rotrV NullV = NullV
-rotrV vs    = tailV vs <: headV vs
-
-rotlV NullV = NullV
-rotlV vs    = lastV vs :> initV vs
-
-rotateV n
-  | n > 0     = pipeV (copyV (abs n) rotlV)
-  | n < 0     = pipeV (copyV (abs n) rotrV)
-  | otherwise = id
-
-
-concatV = foldrV (<+>) NullV
-
-reverseV NullV   = NullV
-reverseV (v:>vs) = reverseV vs <: v
-
-generateV 0 _ _ = NullV
-generateV n f a = x :> generateV (n-1) f x 
-        where x = f a
-
-iterateV 0 _ _ = NullV
-iterateV n f a = a :> iterateV (n-1) f (f a)
-
-copyV k x = iterateV k id x 
-
 {-
 serialV  fs  x = serialV' (reverseV fs ) x
   where
@@ -395,15 +555,6 @@
    = error "parallelV: Vectors have not the same size!"
 parallelV (f:>fs) (x:>xs) = f x :> parallelV fs xs
 
-scanlV _ _ NullV   = NullV
-scanlV f a (x:>xs) = q :> scanlV f q xs 
-       where q = f a x
-
-scanrV _ _ NullV  = NullV
-scanrV f a (x:>NullV) = f x a :> NullV
-scanrV f a (x:>xs)    = f x y :> ys 
-          where ys@(y:>_) = scanrV f a xs
-
 meshlV _ a NullV   = (a, NullV)
 meshlV f a (x:>xs) = (a'', y:>ys) 
        where (a', y)   = f a x
@@ -414,9 +565,4 @@
         where (y, a'') = f x a'
           (ys, a') = meshrV f a xs
 -}
-
-
-
-
-
 
diff --git a/src/ForSyDe/Shallow/MoC/CSDF.hs b/src/ForSyDe/Shallow/MoC/CSDF.hs
--- a/src/ForSyDe/Shallow/MoC/CSDF.hs
+++ b/src/ForSyDe/Shallow/MoC/CSDF.hs
@@ -38,6 +38,9 @@
 -- | The process constructor 'delaynCSDF' delays the signal n event
 --   cycles by introducing n initial values at the beginning of the
 --   output signal.
+--
+-- >>> delayCSDF [3,2,1,0] $ signal [1..5]
+-- {3,2,1,0,1,2,3,4,5}
 delayCSDF :: [a] -> Signal a -> Signal a
 delayCSDF initial_tokens xs = signal initial_tokens +-+ xs
 
@@ -56,6 +59,15 @@
 -- produced tokens and the function) defined in the list of tuples, given as
 -- argument, in a cyclic fashion. The length of the list of scenarios gives the
 -- actor's cycle period.
+--
+-- >>> let c1 = (2,1,\[a,b] -> [a + b])
+-- >>> let c2 = (2,2,\[a,b] -> [max a b, min a b])
+-- >>> let c3 = (2,1,\[a,b] -> [a - b])
+-- >>> let s = signal [1..20]
+-- >>> s
+-- {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}
+-- >>> actor11CSDF [c1,c2,c3] s
+-- {3,4,3,-1,15,10,9,-1,27,16,15,-1,39}
 actor11CSDF :: [(Int, Int, [a] -> [b])] -> Signal a -> Signal b
 actor11CSDF = mapCSDF
 
@@ -65,6 +77,14 @@
 -- produced tokens and the function) defined in the list of tuples, given as
 -- argument, in a cyclic fashion. The length of the list of scenarios gives the
 -- actor's cycle period.
+--
+-- >>> let c1 = ((2,2),1,\[a,b] [c,d]   -> [a+b+c+d])
+-- >>> let c2 = ((1,1),2,\  [a] [b]     -> [max a b, min a b])
+-- >>> let c3 = ((1,3),1,\  [a] [b,c,d] -> [a-b+c-d])
+-- >>> let s1 = signal [1..10]
+-- >>> let s2 = signal [10..] 
+-- >>> actor21CSDF [c1,c2,c3] s1 s2
+-- {24,12,3,-10,44,18,7,-12,64}
 actor21CSDF :: [((Int, Int), Int, [a] -> [b] -> [c])]
             -> Signal a -> Signal b -> Signal c
 actor21CSDF = zipWithCSDF
@@ -107,6 +127,11 @@
 -- produced tokens and the function) defined in the list of tuples, given as
 -- argument, in a cyclic fashion. The length of the list of scenarios gives the
 -- actor's cycle period.
+--
+-- >>> let c1 = ((2,1), (0,1), \[a,b] [c] -> ([], [a+b+c]))
+-- >>> let c2 = ((1,3), (2,3), \[a] [b,c,d] -> ([a,b], [b, c, d]))
+-- >>> actor22CSDF [c1,c2] (signal [1..10]) (signal [11..20])
+-- ({3,12,6,16},{14,12,13,14,24,16,17,18,34})
 actor22CSDF :: [((Int, Int), (Int, Int), [a] -> [b] -> ([c], [d]))]
             -> Signal a -> Signal b -> (Signal c, Signal d)
 actor22CSDF s xs ys = unzipCSDF (outputTokens s) $ zipWithCSDF (inpOut2n s) xs ys
@@ -211,6 +236,11 @@
 -- produced tokens and the function) defined in the list of tuples, given as
 -- argument, in a cyclic fashion. The length of the list of scenarios gives the
 -- actor's cycle period.
+--
+-- >>> let c1 = ((1,0,1), (1,1,3,0), \[a] _ [b] -> ([b], [succ b], [a, 2*a, 3*a], []))
+-- >>> let c2 = ((2,1,1), (0,2,1,1), \[a,b] [c] [d] -> ([], [d, succ d], [a+b], [c]))
+-- >>> actor34CSDF [c1,c2] (signal [1..10]) (signal [11..20]) (signal ['a'..'k'])
+-- ({'a','c','e','g'},{'b','b','c','d','d','e','f','f','g','h'},{1,2,3,5,4,8,12,11,7,14,21,17,10,20,30},{11,12,13})
 actor34CSDF :: [((Int, Int, Int), (Int, Int, Int, Int),
             [a] -> [b] -> [c] -> ([d], [e], [f], [g]))]
             -> Signal a -> Signal b -> Signal c
diff --git a/src/ForSyDe/Shallow/MoC/SADF.hs b/src/ForSyDe/Shallow/MoC/SADF.hs
--- a/src/ForSyDe/Shallow/MoC/SADF.hs
+++ b/src/ForSyDe/Shallow/MoC/SADF.hs
@@ -113,6 +113,13 @@
 -- | The process constructor 'kernel22SADF' constructs an SADF kernel with
 -- two data input and two data output signals. The scenario (token rates and
 -- function) is determined by the control signal.
+--
+-- >>> let scen1 = ((1,1), (1,1), \[a] [b] -> ([2*a], [2*b]))
+-- >>> let scen2 = ((2,2), (1,1), \[a,b] [c,d] -> ([a+b], [c+d]))
+-- >>> let scen3 = ((1,2), (2,1), \[a] [b,c] -> ([b,c], [a]))
+-- >>> let sc = signal [scen1, scen2, scen3]
+-- >>> kernel22SADF sc (signal [1..20]) (signal [21 .. 40])
+-- ({2,5,24,25},{42,45,4}) 
 kernel22SADF :: Signal ((Int, Int), (Int, Int), [a] -> [b] -> ([c], [d]))
              -> Signal a -> Signal b
              -> (Signal c, Signal d)
diff --git a/src/ForSyDe/Shallow/MoC/SDF.hs b/src/ForSyDe/Shallow/MoC/SDF.hs
--- a/src/ForSyDe/Shallow/MoC/SDF.hs
+++ b/src/ForSyDe/Shallow/MoC/SDF.hs
@@ -20,7 +20,7 @@
   -- * Sequential Process Constructors
   -- | Sequential process constructors are used for processes that
   -- have a state. One of the input parameters is the initial state.
-  delaySDF, delaynSDF,
+  delaySDF,
   -- -- * Processes
   -- -- | Processes to unzip a signal of tupels into a tuple of signals
   -- unzipSDF, unzip3SDF, unzip4SDF,
@@ -43,22 +43,16 @@
 -------------------------------------
 
 -- | The process constructor 'delaySDF' delays the signal one event
---   cycle by introducing an initial value at the beginning of the
---   output signal.  Note, that this implies that there is one event
---   (the first) at the output signal that has no corresponding event at
---   the input signal.  One could argue that input and output signals
---   are not fully synchronized, even though all input events are
---   synchronous with a corresponding output event. However, this is
---   necessary to initialize feed-back loops.
-delaySDF :: a -> Signal a -> Signal a
-delaySDF x xs = x :- xs
-
-
--- | The process constructor 'delaynSDF' delays the signal n event
---   cycles by introducing n initial values at the beginning of the
---   output signal.  
-delaynSDF :: [a] -> Signal a -> Signal a 
-delaynSDF initial_tokens xs = signal initial_tokens +-+ xs 
+--   cycle by introducing a set of initial values at the beginning of
+--   the output signal.  Note, that this implies that there is a
+--   prefix at the output signal (the first n events) that has no
+--   corresponding event at the input signal. This is necessary to
+--   initialize feedback loops.
+--
+-- >>> delaySDF [0,0,0] $ signal [1,2,3,4]
+-- {0,0,0,1,2,3,4}
+delaySDF :: [a] -> Signal a -> Signal a 
+delaySDF initial_tokens xs = signal initial_tokens +-+ xs 
 
 ------------------------------------------------------------------------
 --
@@ -72,15 +66,25 @@
 -- one input and one output signals. For each input or output signal,
 -- the process constructor takes the number of consumed and produced
 -- tokens and the function of the actor as arguments.
+--
+-- >>> let f [a,b] = [a+b,a-b,a*b]
+-- >>> actor11SDF 2 3 f $ signal [1,2,3,4,5]
+-- {3,-1,2,7,-1,12}
 actor11SDF :: Int -> Int -> ([a] -> [b]) -> Signal a -> Signal b
-actor11SDF c p f s1 = produceSDF p $ mapSDF c f s1 
+actor11SDF = mapSDF 
 
 -- | The process constructor 'actor21SDF' constructs an SDF actor with
 -- two input and one output signals. For each input or output signal,
 -- the process constructor takes the number of consumed and produced
 -- tokens and the function of the actor as arguments.
+--
+-- >>> let f [a,b] [c] = [a+b+c,b-c]
+-- >>> let s1 = signal [1..6]
+-- >>> let s2 = signal [1..]
+-- >>> actor21SDF (2,1) 2 f s1 s2
+-- {4,1,9,2,14,3}
 actor21SDF :: (Int, Int) -> Int -> ([a] -> [b] -> [c]) -> Signal a -> Signal b -> Signal c    
-actor21SDF c p f s1 s2 = produceSDF p $ zipWithSDF c f s1 s2
+actor21SDF = zipWithSDF
 
 -- | The process constructor 'actor31SDF' constructs an SDF actor with
 -- three input and one output signals. For each input or output signal,
@@ -88,7 +92,7 @@
 -- tokens and the function of the actor as arguments.
 actor31SDF :: (Int, Int, Int) -> Int -> ([a] -> [b] -> [c] -> [d])
        -> Signal a -> Signal b -> Signal c -> Signal d   
-actor31SDF c p f s1 s2 s3 = produceSDF p $ zipWith3SDF c f s1 s2 s3
+actor31SDF = zipWith3SDF
 
 -- | The process constructor 'actor41SDF' constructs an SDF actor with
 -- four input and one output signals. For each input or output signal,
@@ -97,7 +101,7 @@
 actor41SDF :: (Int, Int, Int, Int) -> Int 
     -> ([a] -> [b] -> [c] -> [d] -> [e]) 
     -> Signal a -> Signal b -> Signal c -> Signal d -> Signal e 
-actor41SDF c p f s1 s2 s3 s4 = produceSDF p $ zipWith4SDF c f s1 s2 s3 s4
+actor41SDF = zipWith4SDF
 
 
 -- > Actors with two outputs
@@ -108,7 +112,7 @@
 -- tokens and the function of the actor as arguments.
 actor12SDF :: Int -> (Int, Int) -> ([a] -> ([b], [c]))
            -> Signal a -> (Signal b, Signal c)
-actor12SDF c (p1,p2) f xs = unzipSDF (p1, p2) $ mapSDF c f xs  
+actor12SDF c (p1,p2) f xs = unzipSDF (p1,p2) $ mapSDF c 1 (wrapF1 f) xs  
 
 -- | The process constructor 'actor22SDF' constructs an SDF actor with
 -- two input and two output signals. For each input or output signal,
@@ -116,7 +120,7 @@
 -- tokens and the function of the actor as arguments.
 actor22SDF :: (Int, Int) -> (Int, Int) -> ([a] -> [b] -> ([c], [d]))
            -> Signal a -> Signal b -> (Signal c, Signal d)
-actor22SDF (c1, c2) (p1, p2) f xs ys = unzipSDF (p1, p2) $ zipWithSDF (c1, c2) f xs ys
+actor22SDF (c1,c2) (p1,p2) f xs ys = unzipSDF (p1,p2) $ zipWithSDF (c1,c2) 1 (wrapF2 f) xs ys
 
 -- | The process constructor 'actor32SDF' constructs an SDF actor with
 -- three input and two output signals. For each input or output signal,
@@ -125,8 +129,8 @@
 actor32SDF :: (Int, Int, Int) -> (Int, Int)
            -> ([a] -> [b] -> [c] -> ([d], [e]))
            -> Signal a -> Signal b -> Signal c -> (Signal d, Signal e)
-actor32SDF (c1, c2, c3) (p1, p2) f as bs cs
-  = unzipSDF (p1, p2) $ zipWith3SDF (c1, c2, c3) f as bs cs
+actor32SDF (c1,c2,c3) (p1,p2) f as bs cs
+  = unzipSDF (p1,p2) $ zipWith3SDF (c1,c2,c3) 1 (wrapF3 f) as bs cs
 
 -- | The process constructor 'actor42SDF' constructs an SDF actor with
 -- four input and two output signals. For each input or output signal,
@@ -136,8 +140,8 @@
            -> ([a] -> [b] -> [c] -> [d] -> ([e], [f])) 
            -> Signal a -> Signal b -> Signal c -> Signal d 
            -> (Signal e, Signal f)
-actor42SDF (c1, c2, c3, c4) (p1, p2) f as bs cs ds 
-  = unzipSDF (p1, p2)$ zipWith4SDF (c1, c2, c3, c4) f as bs cs ds
+actor42SDF (c1,c2,c3,c4) (p1,p2) f as bs cs ds 
+  = unzipSDF (p1,p2) $ zipWith4SDF (c1,c2,c3,c4) 1 (wrapF4 f) as bs cs ds
 
 -- > Actors with three outputs
 
@@ -148,7 +152,7 @@
 actor13SDF :: Int -> (Int, Int, Int) 
            -> ([a] -> ([b], [c], [d])) 
            -> Signal a -> (Signal b, Signal c, Signal d)
-actor13SDF c (p1, p2, p3) f xs = unzip3SDF (p1, p2, p3) $ mapSDF c f xs  
+actor13SDF c (p1,p2,p3) f xs = unzip3SDF (p1,p2,p3) $ mapSDF c 1 (wrapF1 f) xs  
 
 -- | The process constructor 'actor23SDF' constructs an SDF actor with
 -- two input and three output signals. For each input or output signal,
@@ -158,8 +162,8 @@
            -> ([a] -> [b] -> ([c], [d], [e])) 
            -> Signal a -> Signal b 
            -> (Signal c, Signal d, Signal e)
-actor23SDF (c1, c2) (p1, p2, p3) f xs ys
-  = unzip3SDF (p1, p2, p3) $ zipWithSDF (c1, c2) f xs ys
+actor23SDF (c1,c2) (p1,p2,p3) f xs ys
+  = unzip3SDF (p1,p2,p3) $ zipWithSDF (c1,c2) 1 (wrapF2 f) xs ys
 
 -- | The process constructor 'actor33SDF' constructs an SDF actor with
 -- three input and three output signals. For each input or output signal,
@@ -168,8 +172,8 @@
 actor33SDF :: (Int, Int, Int) -> (Int, Int, Int) 
            -> ([a] -> [b] -> [c] -> ([d], [e], [f])) 
            -> Signal a -> Signal b -> Signal c -> (Signal d, Signal e, Signal f)
-actor33SDF (c1, c2, c3) (p1, p2, p3) f as bs cs
-  = unzip3SDF (p1, p2, p3) $ zipWith3SDF (c1, c2, c3) f as bs cs
+actor33SDF (c1,c2,c3) (p1,p2,p3) f as bs cs
+  = unzip3SDF (p1,p2,p3) $ zipWith3SDF (c1,c2,c3) 1 (wrapF3 f) as bs cs
 
 -- | The process constructor 'actor43SDF' constructs an SDF actor with
 -- four input and three output signals. For each input or output signal,
@@ -179,8 +183,8 @@
            -> ([a] -> [b] -> [c] -> [d] -> ([e], [f], [g])) 
            -> Signal a -> Signal b -> Signal c -> Signal d 
            -> (Signal e, Signal f, Signal g)
-actor43SDF (c1, c2, c3, c4) (p1, p2, p3) f as bs cs ds 
-  = unzip3SDF (p1, p2, p3)$ zipWith4SDF (c1, c2, c3, c4) f as bs cs ds
+actor43SDF (c1,c2,c3,c4) (p1,p2,p3) f as bs cs ds 
+  = unzip3SDF (p1,p2,p3)$ zipWith4SDF (c1,c2,c3,c4) 1 (wrapF4 f) as bs cs ds
 
 -- > Actors with four outputs
 
@@ -191,7 +195,7 @@
 actor14SDF :: Int -> (Int, Int, Int, Int) 
            -> ([a] -> ([b], [c], [d], [e])) 
            -> Signal a -> (Signal b, Signal c, Signal d, Signal e)
-actor14SDF c (p1, p2, p3, p4) f xs = unzip4SDF (p1, p2, p3, p4) $ mapSDF c f xs  
+actor14SDF c (p1,p2,p3,p4) f xs = unzip4SDF (p1,p2,p3,p4) $ mapSDF c 1 (wrapF1 f) xs  
 
 -- | The process constructor 'actor24SDF' constructs an SDF actor with
 -- two input and four output signals. For each input or output signal,
@@ -201,8 +205,8 @@
        -> ([a] -> [b] -> ([c], [d], [e], [f])) 
        -> Signal a -> Signal b 
        -> (Signal c, Signal d, Signal e, Signal f)
-actor24SDF (c1, c2) (p1, p2, p3, p4) f xs ys
-  = unzip4SDF (p1, p2, p3, p4) $ zipWithSDF (c1, c2) f xs ys
+actor24SDF (c1,c2) (p1,p2,p3,p4) f xs ys
+  = unzip4SDF (p1,p2,p3,p4) $ zipWithSDF (c1,c2) 1 (wrapF2 f) xs ys
 
 -- | The process constructor 'actor34SDF' constructs an SDF actor with
 -- three input and four output signals. For each input or output signal,
@@ -212,8 +216,8 @@
            -> ([a] -> [b] -> [c] -> ([d], [e], [f], [g])) 
            -> Signal a -> Signal b -> Signal c
            -> (Signal d, Signal e, Signal f, Signal g)
-actor34SDF (c1, c2, c3) (p1, p2, p3, p4) f as bs cs 
-  = unzip4SDF (p1, p2, p3, p4) $ zipWith3SDF (c1, c2, c3) f as bs cs
+actor34SDF (c1,c2,c3) (p1,p2,p3,p4) f as bs cs 
+  = unzip4SDF (p1,p2,p3,p4) $ zipWith3SDF (c1,c2,c3) 1 (wrapF3 f) as bs cs
 
 -- | The process constructor 'actor14SDF' constructs an SDF actor with
 -- four input and four output signals. For each input or output signal,
@@ -223,8 +227,8 @@
            -> ([a] -> [b] -> [c] -> [d] -> ([e], [f], [g], [h])) 
            -> Signal a -> Signal b -> Signal c -> Signal d 
            -> (Signal e, Signal f, Signal g, Signal h)
-actor44SDF (c1, c2, c3, c4) (p1, p2, p3, p4) f as bs cs ds 
-  = unzip4SDF (p1, p2, p3, p4)$ zipWith4SDF (c1, c2, c3, c4) f as bs cs ds
+actor44SDF (c1,c2,c3,c4) (p1,p2,p3,p4) f as bs cs ds 
+  = unzip4SDF (p1,p2,p3,p4) $ zipWith4SDF (c1,c2,c3,c4) 1 (wrapF4 f) as bs cs ds
 
 ---------------------------------------------------------------------
 -- COMBINATIONAL PROCESS CONSTRUCTORS
@@ -234,12 +238,15 @@
 -- (@c@) and produced (@p@) tokens and a function @f@ that operates on
 -- a list, and results in an SDF-process that takes an input signal
 -- and results in an output signal
-mapSDF :: Int -> ([a] -> b) -> Signal a -> Signal b
-mapSDF _ _ NullS   = NullS
-mapSDF c f xs     
+mapSDF :: Int -> Int -> ([a] -> [b]) -> Signal a -> Signal b
+mapSDF _ _ _ NullS   = NullS
+mapSDF c p f xs     
   | c <= 0 = error "mapSDF: Number of consumed tokens must be positive integer" 
   | not $ sufficient_tokens c xs  = NullS
-  | otherwise  = produced_tokens :- mapSDF c f (dropS c xs) 
+  | otherwise  = if length produced_tokens == p then
+                   signal produced_tokens +-+ mapSDF c p f (dropS c xs) 
+                 else   
+                   error "mapSDF: Function does not produce correct number of tokens" 
   where consumed_tokens = fromSignal $ takeS c xs
         produced_tokens = f consumed_tokens
 
@@ -248,15 +255,18 @@
 -- the number of produced tokens and a function @f@
 -- that operates on two lists, and results in an SDF-process that takes two
 -- input signals and results in an output signal
-zipWithSDF :: (Int, Int) -> ([a] -> [b] -> c)
+zipWithSDF :: (Int, Int) -> Int -> ([a] -> [b] -> [c])
            -> Signal a -> Signal b -> Signal c                  
-zipWithSDF (_, _) _ NullS _ = NullS
-zipWithSDF (_, _) _ _ NullS = NullS
-zipWithSDF (c1, c2) f as bs 
+zipWithSDF (_, _) _ _ NullS _ = NullS
+zipWithSDF (_, _) _ _ _ NullS = NullS
+zipWithSDF (c1, c2) p f as bs 
   | c1 <= 0 || c2 <= 0  = error "zipWithSDF: Number of consumed tokens must be positive integer"
   | (not $ sufficient_tokens c1 as) 
     || (not $ sufficient_tokens c2 bs) = NullS
-  | otherwise = produced_tokens :- zipWithSDF (c1, c2) f (dropS c1 as) (dropS c2 bs)  
+  | otherwise = if length produced_tokens == p then
+                  signal produced_tokens +-+ zipWithSDF (c1, c2) p f (dropS c1 as) (dropS c2 bs)  
+                else
+                  error "zipWithSDF: Function does not produce correct number of tokens"
   where consumed_tokens_as = fromSignal $ takeS c1 as
         consumed_tokens_bs = fromSignal $ takeS c2 bs
         produced_tokens = f consumed_tokens_as consumed_tokens_bs
@@ -266,12 +276,12 @@
 -- the number of produced tokens and a function @f@
 -- that operates on three lists, and results in an SDF-process that takes three
 -- input signals and results in an output signal  
-zipWith3SDF :: (Int, Int, Int) -> ([a] -> [b] -> [c] -> d) 
+zipWith3SDF :: (Int, Int, Int) -> Int -> ([a] -> [b] -> [c] -> [d]) 
             -> Signal a -> Signal b -> Signal c -> Signal d                 
-zipWith3SDF (_, _, _) _ NullS _ _= NullS
-zipWith3SDF (_, _, _) _ _ NullS _= NullS
-zipWith3SDF (_, _, _) _ _ _ NullS= NullS
-zipWith3SDF (c1, c2, c3) f as bs cs
+zipWith3SDF (_, _, _) _ _ NullS _ _= NullS
+zipWith3SDF (_, _, _) _ _ _ NullS _= NullS
+zipWith3SDF (_, _, _) _ _ _ _ NullS= NullS
+zipWith3SDF (c1, c2, c3) p f as bs cs
   | c1 <= 0 || c2 <= 0 || c3 <= 0
   = error "zipWith3SDF: Number of consumed tokens must be positive integer"
   | (not $ sufficient_tokens c1 as) 
@@ -279,8 +289,11 @@
     || (not $ sufficient_tokens c3 cs)
   = NullS
   | otherwise
-  = produced_tokens :- zipWith3SDF (c1, c2, c3) f
-                       (dropS c1 as) (dropS c2 bs) (dropS c3 cs)
+  = if length produced_tokens == p then
+      signal produced_tokens +-+ zipWith3SDF (c1, c2, c3) p f
+                                 (dropS c1 as) (dropS c2 bs) (dropS c3 cs)
+    else
+      error "zipWith3SDF: Function does not produce correct number of tokens"
   where consumed_tokens_as = fromSignal $ takeS c1 as
         consumed_tokens_bs = fromSignal $ takeS c2 bs
         consumed_tokens_cs = fromSignal $ takeS c3 cs
@@ -292,14 +305,14 @@
 -- denoting the number of produced tokens and a function @f@ that
 -- operates on three lists, and results in an SDF-process that takes
 -- three input signals and results in an output signal
-zipWith4SDF :: (Int, Int, Int, Int) 
-            -> ([a] -> [b] -> [c] -> [d] -> e) 
+zipWith4SDF :: (Int, Int, Int, Int) -> Int 
+            -> ([a] -> [b] -> [c] -> [d] -> [e]) 
             -> Signal a -> Signal b -> Signal c -> Signal d -> Signal e 
-zipWith4SDF (_, _, _, _) _ NullS _ _ _ = NullS
-zipWith4SDF (_, _, _, _) _ _ NullS _ _ = NullS
-zipWith4SDF (_, _, _, _) _ _ _ NullS _ = NullS
-zipWith4SDF (_, _, _, _) _ _ _ _ NullS = NullS
-zipWith4SDF (c1, c2, c3, c4) f as bs cs ds
+zipWith4SDF (_, _, _, _) _ _ NullS _ _ _ = NullS
+zipWith4SDF (_, _, _, _) _ _ _ NullS _ _ = NullS
+zipWith4SDF (_, _, _, _) _ _ _ _ NullS _ = NullS
+zipWith4SDF (_, _, _, _) _ _ _ _ _ NullS = NullS
+zipWith4SDF (c1, c2, c3, c4) p f as bs cs ds
   | c1 <= 0 || c2 <= 0 || c3 <= 0 || c4 <= 0
   = error "zipWith4SDF: Number of consumed tokens must be positive integer"
   | (not $ sufficient_tokens c1 as) 
@@ -308,26 +321,22 @@
     || (not $ sufficient_tokens c4 ds)    
   = NullS
   | otherwise    
-  = produced_tokens :- zipWith4SDF (c1, c2, c3, c4) f
-           (dropS c1 as) (dropS c2 bs) (dropS c3 cs) (dropS c4 ds)
+  = if length produced_tokens == p then
+      signal produced_tokens +-+ zipWith4SDF (c1, c2, c3, c4) p f
+             (dropS c1 as) (dropS c2 bs) (dropS c3 cs) (dropS c4 ds)
+    else
+      error "zipWith4SDF: Function does not produce correct number of tokens"
   where consumed_tokens_as = fromSignal $ takeS c1 as
         consumed_tokens_bs = fromSignal $ takeS c2 bs
         consumed_tokens_cs = fromSignal $ takeS c3 cs
         consumed_tokens_ds = fromSignal $ takeS c4 ds
         produced_tokens = f consumed_tokens_as consumed_tokens_bs
                             consumed_tokens_cs consumed_tokens_ds
-        
-
+                            
 ---------------------------------------------------------------------
 -- unzipSDF Processes
 ---------------------------------------------------------------------
 
-produceSDF :: Int -> Signal [a] -> Signal a
-produceSDF p s
-  | allS (\partition -> length partition == p) s
-  = foldrS (\part sig -> signal part +-+ sig) NullS s
-  | otherwise = error "SDF: Function does not produce correct number of tokens"
-
 unzipSDF :: (Int, Int) -> Signal ([a], [b]) 
          -> (Signal a, Signal b)
 unzipSDF (p1, p2) xs = (s1, s2) 
@@ -336,16 +345,12 @@
     s2 = signal $ f2 xs
     f1 NullS     = []
     f1 ((as, _):-xs)
-      = if length as == p1 then 
-          as ++ f1 xs
-        else  
-          error "unzipSDF: Process does not produce correct number of tokens"
+      | length as == p1 = as ++ f1 xs
+      | otherwise = error "unzipSDF: Process does not produce correct number of tokens"
     f2 NullS     = []
     f2 ((_, bs):-xs)
-      = if length bs == p2 then 
-          bs ++ f2 xs
-        else  
-          error "unzipSDF: Process does not produce correct number of tokens"  
+      | length bs == p2 = bs ++ f2 xs
+      | otherwise = error "unzipSDF: Process does not produce correct number of tokens"  
 
 
 unzip3SDF :: (Int, Int, Int) -> Signal ([a], [b], [c]) 
@@ -357,22 +362,16 @@
     s3 = signal $ f3 xs
     f1 NullS      = []
     f1 ((as, _, _):-xs) 
-      = if length as == p1 then
-          as ++ f1 xs
-        else  
-          error "unzip3SDF: Process does not produce correct number of tokens"
+      | length as == p1 = as ++ f1 xs
+      | otherwise = error "unzip3SDF: Process does not produce correct number of tokens"
     f2 NullS      = []
     f2 ((_, bs, _):-xs) 
-      = if length bs == p2 then 
-          bs ++ f2 xs
-        else  
-          error "unzip3SDF: Process does not produce correct number of tokens"  
+      | length bs == p2 = bs ++ f2 xs
+      | otherwise = error "unzip3SDF: Process does not produce correct number of tokens"  
     f3 NullS      = []
     f3 ((_, _, cs):-xs) 
-      = if length cs == p3 then 
-          cs ++ f3 xs
-        else  
-          error "unzip3SDF: Process does not produce correct number of tokens"  
+      | length cs == p3 = cs ++ f3 xs
+      | otherwise = error "unzip3SDF: Process does not produce correct number of tokens"  
   
 
 unzip4SDF :: (Int, Int, Int, Int) -> Signal ([a], [b], [c], [d]) 
@@ -385,28 +384,20 @@
     s4 = signal $ f4 xs
     f1 NullS      = []
     f1 ((as, _, _, _):-xs) 
-      = if length as == p1 then
-          as ++ f1 xs
-        else  
-          error "unzip4SDF: Process does not produce correct number of tokens"
+      | length as == p1 = as ++ f1 xs
+      | otherwise = error "unzip4SDF: Process does not produce correct number of tokens"
     f2 NullS      = []
     f2 ((_, bs, _, _):-xs) 
-      = if length bs == p2 then 
-          bs ++ f2 xs
-        else  
-          error "unzip4SDF: Process does not produce correct number of tokens"  
+      | length bs == p2 = bs ++ f2 xs
+      | otherwise = error "unzip4SDF: Process does not produce correct number of tokens"  
     f3 NullS      = []
     f3 ((_, _, cs, _):-xs) 
-      = if length cs == p3 then 
-          cs ++ f3 xs
-        else  
-          error "unzip4SDF: Process does not produce correct number of tokens" 
+      | length cs == p3 = cs ++ f3 xs
+      | otherwise = error "unzip4SDF: Process does not produce correct number of tokens" 
     f4 NullS      = []
     f4 ((_, _, _, ds):-xs) 
-      = if length ds == p4 then 
-          ds ++ f4 xs
-        else  
-          error "unzip4SDF: Process does not produce correct number of tokens" 
+      | length ds == p4 = ds ++ f4 xs
+      | otherwise = error "unzip4SDF: Process does not produce correct number of tokens" 
 
 ------------------------------------------------------------------------
 --
@@ -414,14 +405,31 @@
 --
 ------------------------------------------------------------------------
 
+
 sufficient_tokens :: (Num a, Eq a, Ord a) => a -> Signal t -> Bool
 sufficient_tokens 0 _     = True
 sufficient_tokens _ NullS = False
 sufficient_tokens n (_:-xs)
-  = if n < 0 then
-      error "sufficient_tokens: n must not be negative"
-    else
-      sufficient_tokens (n-1) xs
+ = if n < 0 then
+     error "sufficient_tokens: n must not be negative"
+   else
+     sufficient_tokens (n-1) xs
+
+-- outputTokens :: [(a, b, c)] -> [b]
+-- outputTokens [] = []
+-- outputTokens ((_, b, _):xs) = b : outputTokens xs
+
+wrapF1 :: ([a] -> y) -> [a] -> [y]
+wrapF1 f = (\a -> [f a])
+
+wrapF2 :: ([a] -> [b] -> y) -> [a] -> [b] -> [y]
+wrapF2 f = (\a b -> [f a b])
+
+wrapF3 :: ([a] -> [b] -> [c] -> y) -> [a] -> [b] -> [c] -> [y]
+wrapF3 f = (\a b c -> [f a b c])
+
+wrapF4 :: ([a] -> [b] -> [c] -> [d] -> y) -> [a] -> [b] -> [c] -> [d] -> [y]
+wrapF4 f = (\a b c d -> [f a b c d])
 
 ------------------------------------------------------------------------
 --
diff --git a/src/ForSyDe/Shallow/MoC/Synchronous/Lib.hs b/src/ForSyDe/Shallow/MoC/Synchronous/Lib.hs
--- a/src/ForSyDe/Shallow/MoC/Synchronous/Lib.hs
+++ b/src/ForSyDe/Shallow/MoC/Synchronous/Lib.hs
@@ -16,7 +16,7 @@
   -- | Combinational process constructors are used for processes that
   --   do not have a state.
   mapSY, zipWithSY, zipWith3SY, 
-  zipWith4SY, zipWithxSY,
+  zipWith4SY, mapxSY, zipWithxSY,
   combSY, comb2SY, comb3SY, comb4SY,
   -- ** Sequential process constructors
   -- | Sequential process constructors are used for processes that
@@ -31,7 +31,7 @@
   --   to many cases.
   whenSY, zipSY, zip3SY, zip4SY, zip5SY, zip6SY, 
   unzipSY, unzip3SY, unzip4SY, unzip5SY, unzip6SY,
-  zipxSY, unzipxSY, mapxSY, 
+  zipxSY, unzipxSY, 
   fstSY, sndSY
   ) where
 
@@ -44,7 +44,9 @@
 -- | The process constructor 'mapSY' takes a combinational function as
 -- argument and returns a process with one input signal and one output
 -- signal.
-
+--
+-- >>> mapSY (+1) $ signal [1,2,3,4]
+-- {2,3,4,5}
 mapSY :: (a -> b) -> Signal a -> Signal b
 mapSY _ NullS   = NullS
 mapSY f (x:-xs) = f x :- (mapSY f xs)
@@ -52,6 +54,9 @@
 -- | The process constructor 'zipWithSY' takes a combinational
 -- function as argument and returns a process with two input signals
 -- and one output signal.
+--
+-- >>> zipWithSY (+) (signal [1,2,3,4]) (signal [11,12,13,14,15,16,17])
+-- {12,14,16,18}
 zipWithSY :: (a -> b -> c) -> Signal a -> Signal b -> Signal c
 zipWithSY _ NullS   _   = NullS
 zipWithSY _ _   NullS   = NullS
@@ -101,12 +106,24 @@
 comb4SY = zipWith4SY
    
 -- | The process constructor 'mapxSY' creates a process network that
--- maps a function onto all signals in a vector of signals.
+-- maps a function onto all signals in a vector of signals. See 'mapV'.
+--
+-- >>> let s1 = signal [1,2,3,4]
+-- >>> let s2 = signal [10,20,30,40]
+-- >>> let s3 = signal [100,200,300]
+-- >>> mapxSY (+1) $ vector [s1,s2,s3]
+-- <{2,3,4,5},{11,21,31,41},{101,201,301}> 
 mapxSY :: (a -> b) -> Vector (Signal a) -> Vector (Signal b)
 mapxSY f = mapV (mapSY f)
 
 -- | The process constructor 'zipWithxSY' works as 'zipWithSY', but
 -- takes a vector of signals as input.
+--
+-- >>> let s1 = signal [1,2,3,4]
+-- >>> let s2 = signal [10,20,30,40]
+-- >>> let s3 = signal [100,200,300]
+-- >>> zipWithxSY (reduceV (+)) $ vector [s1,s2,s3]
+-- {111,222,333}
 zipWithxSY :: (Vector a -> b) -> Vector (Signal a) -> Signal b
 zipWithxSY f = mapSY f . zipxSY
 
@@ -122,6 +139,9 @@
 -- are not fully synchronized, even though all input events are
 -- synchronous with a corresponding output event. However, this is
 -- necessary to initialize feed-back loops.
+--
+-- >>> delaySY 1 $ signal [1,2,3,4]
+-- {1,1,2,3,4}
 delaySY :: a        -- ^Initial state
         -> Signal a -- ^Input signal
         -> Signal a -- ^Output signal
@@ -129,6 +149,9 @@
 
 -- | The process constructor 'delaynSY' delays the signal n events by
 -- introducing n identical default values.
+--
+-- >>> delaynSY 0 3 $ signal [1,2,3,4]
+-- {0,0,0,1,2,3,4}
 delaynSY :: a        -- ^Initial state
          -> Int      -- ^ Delay cycles 
          -> Signal a -- ^Input signal
@@ -143,9 +166,8 @@
 -- 'scanlSY' and has the value of the new state as its output value as
 -- illustrated by the following example.
 --
--- > SynchronousLib> scanlSY (+) 0 (signal [1,2,3,4])
---
--- > {1,3,6,10} :: Signal Integer
+-- >>> scanlSY (+) 0 (signal [1,2,3,4])
+-- {1,3,6,10}
 -- 
 -- This is in contrast to the function 'scanldSY', which has its
 -- current state as its output value.
@@ -178,9 +200,8 @@
 -- the output value is the current state and not the one of the next
 -- state.
 --
--- > SynchronousLib> scanldSY (+) 0 (signal [1,2,3,4])
---
--- > {0,1,3,6,10} :: Signal Integer
+-- >>> scanldSY (+) 0 (signal [1,2,3,4])
+-- {0,1,3,6,10}
 scanldSY :: (a -> b -> a) -- ^Combinational function for next state
                           -- decoder
     -> a                  -- ^Initial state
@@ -212,6 +233,9 @@
 -- for the initial state.
 --
 -- In contrast the output of a process created by the process constructor 'mealySY' depends not only on the state, but also on the input values.
+--
+-- >>> mooreSY (+) (*2) 0 $ signal [1,2,3,4,5]
+-- {0,2,6,12,20,30}
 mooreSY :: (a -> b -> a) -- ^Combinational function for next state
                          -- decoder
         -> (a -> c)      -- ^Combinational function for output decoder
@@ -247,6 +271,9 @@
 -- In contrast the output of a process created by the process
 -- constructor 'mooreSY' depends only on the state, but not on the
 -- input values.
+--
+-- >>> mealySY (+) (+) 0 $ signal [1,2,3,4,5]
+-- {1,3,6,10,15}
 mealySY :: (a -> b -> a) -- ^Combinational function for next state
                          -- decoder
        -> (a -> b -> c)  -- ^Combinational function for output decoder
@@ -291,9 +318,8 @@
 -- The process that has the infinite signal of natural numbers as
 -- output is constructed by
 --
--- > SynchronousLib> takeS 5 (sourceSY (+1) 0)
---
--- > {0,1,2,3,4} :: Signal Integer
+-- >>> takeS 5 $ sourceSY (+1) 0
+-- {0,1,2,3,4}
 sourceSY :: (a -> a) -> a -> Signal a
 sourceSY f s0 = o
    where o = delaySY s0 s
@@ -302,6 +328,10 @@
 -- | The process constructor 'fillSY' creates a process that 'fills' a
 -- signal with present values by replacing absent values with a given
 -- value. The output signal is not any more of the type 'AbstExt'.
+--
+-- >>> let s = signal [Abst, Prst 1, Prst 2, Abst, Abst, Prst 4, Abst]
+-- >>> fillSY 3 s
+-- {3,1,2,3,3,4,3}
 fillSY :: a                  -- ^Default value
        -> Signal (AbstExt a) -- ^Absent extended input signal
        -> Signal a           -- ^Output signal
@@ -310,6 +340,10 @@
             replaceAbst _  (Prst x) = x
 
 -- | The process constructor 'holdSY' creates a process that 'fills' a signal with values by replacing absent values by the preceding present value. Only in cases, where no preceding value exists, the absent value is replaced by a default value. The output signal is not any more of the type 'AbstExt'.
+--
+-- >>> let s = signal [Abst, Prst 1, Prst 2, Abst, Abst, Prst 4, Abst]
+-- >>> holdSY 3 s
+-- {3,1,2,2,2,4,4}
 holdSY :: a                  -- ^Default value
        -> Signal (AbstExt a) -- ^Absent extended input signal
        -> Signal a           -- ^Output signal
@@ -326,6 +360,11 @@
 -- of absent extended values. The output signal has the value of the
 -- first signal whenever an event has a present value and 'Abst' when
 -- the event has an absent value.
+--
+-- >>> let clk = signal [Abst,    Prst (), Prst (), Abst,    Abst,    Prst (), Abst]
+-- >>> let sig = signal [Prst 10, Prst 11, Prst 12, Prst 13, Prst 14, Prst 15]
+-- >>> sig `whenSY` clk
+-- {_,11,12,_,_,15}
 whenSY :: Signal (AbstExt a) -> Signal (AbstExt b) 
        -> Signal (AbstExt a)
 whenSY NullS   _          = NullS
@@ -335,12 +374,16 @@
 
 -- | The process 'zipSY' \"zips\" two incoming signals into one signal
 -- of tuples.
+--
+-- >>> zipSY (signal [1,2,3,4]) (signal [10,11,12,13,14])
+-- {(1,10),(2,11),(3,12),(4,13)}
 zipSY :: Signal a -> Signal b -> Signal (a,b)
 zipSY (x:-xs) (y:-ys) = (x, y) :- zipSY xs ys
 zipSY _       _       = NullS
 
 -- | The process 'zip3SY' works as 'zipSY', but takes three input
 -- signals.
+-- 
 zip3SY :: Signal a -> Signal b -> Signal c -> Signal (a,b,c)
 zip3SY (x:-xs) (y:-ys) (z:-zs) = (x, y, z) :- zip3SY xs ys zs
 zip3SY _       _       _       = NullS
@@ -408,13 +451,23 @@
   = (x1:-x1s, x2:-x2s, x3:-x3s, x4:-x4s, x5:-x5s, x6:-x6s)
   where (x1s, x2s, x3s, x4s, x5s, x6s) = unzip6SY xs
 
--- | The process 'zipxSY' \"zips\" a signal of vectors into a vector
--- of signals.
+-- | The process 'zipxSY' transposes a signal of vectors into a vector
+-- of signals. All the events carried by the output signal are synchronized values from all input signals.
+--
+-- >>> let s1 = signal [1,2,3,4]
+-- >>> let s2 = signal [10,20,30,40]
+-- >>> let s3 = signal [100,200,300]
+-- >>> zipxSY $ vector [s1,s2,s3]
+-- {<1,10,100>,<2,20,200>,<3,30,300>}
 zipxSY :: Vector (Signal a) -> Signal (Vector a)
-zipxSY NullV            = NullS
-zipxSY (NullS   :> xss) = zipxSY xss
-zipxSY ((x:-xs) :> xss) = (x :> (mapV headS xss)) 
-                          :- (zipxSY (xs :> (mapV tailS xss)))
+zipxSY = reduceV (zipWithSY (<+>)) . mapV (mapSY unitV)
+  
+-- unsafe implementation! 
+-- zipxSY :: Vector (Signal a) -> Signal (Vector a)
+-- zipxSY NullV            = NullS
+-- zipxSY (NullS   :> xss) = zipxSY xss
+-- zipxSY ((x:-xs) :> xss) = (x :> (mapV headS xss)) 
+--                           :- (zipxSY (xs :> (mapV tailS xss)))
 
 -- | The process 'unzipxSY' \"unzips\" a vector of signals into a
 -- signal of vectors.
diff --git a/src/ForSyDe/Shallow/Utility.hs b/src/ForSyDe/Shallow/Utility.hs
--- a/src/ForSyDe/Shallow/Utility.hs
+++ b/src/ForSyDe/Shallow/Utility.hs
@@ -16,7 +16,9 @@
 module ForSyDe.Shallow.Utility (  
   module ForSyDe.Shallow.Utility.DFT,    
   module ForSyDe.Shallow.Utility.Memory,
-  module ForSyDe.Shallow.Utility.Queue
+  module ForSyDe.Shallow.Utility.Queue,
+  module ForSyDe.Shallow.Utility.Matrix
+  --module ForSyDe.Shallow.Utility.BitVector,
   --module ForSyDe.Shallow.Utility.FilterLib,
   --module ForSyDe.Shallow.Utility.Gaussian,
   --module ForSyDe.Shallow.Utility.PolyArith,
@@ -26,6 +28,8 @@
 import ForSyDe.Shallow.Utility.DFT      
 import ForSyDe.Shallow.Utility.Memory
 import ForSyDe.Shallow.Utility.Queue
+import ForSyDe.Shallow.Utility.Matrix
+--import ForSyDe.Shallow.Utility.BitVector
 --import ForSyDe.Shallow.Utility.FilterLib
 --import ForSyDe.Shallow.Utility.Gaussian
 --import ForSyDe.Shallow.Utility.PolyArith
diff --git a/src/ForSyDe/Shallow/Utility/BitVector.hs b/src/ForSyDe/Shallow/Utility/BitVector.hs
new file mode 100644
--- /dev/null
+++ b/src/ForSyDe/Shallow/Utility/BitVector.hs
@@ -0,0 +1,122 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module  :  ForSyDe.Shallow.Utility.BitVector
+-- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
+-- License     :  BSD-style (see the file LICENSE)
+-- 
+-- Maintainer  :  forsyde-dev@ict.kth.se
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- It defines the bit vector operations from\/to integer.
+-----------------------------------------------------------------------------
+module ForSyDe.Shallow.Utility.BitVector(
+    -- *Polynomial data type
+    BitVector, Parity(..),
+    -- *Bit-vector and integer transformations
+    intToBitVector, bitVectorToInt,
+    -- *Add even\/odd parity bit
+    addEvenParityBit, addOddParityBit, addParityBit,
+    -- *Remove parity bit
+    removeParityBit,
+    -- *Check the even\/odd parity of the bit-vector
+    isEvenParity, isOddParity,
+    -- *Judge the form sanity of the bit-vector
+    isBitVector
+  )
+  where
+
+import ForSyDe.Shallow.Core.Vector
+
+type BitVector = Vector Integer
+
+-- |To judge whether the input bit-vector is in a proper form.
+isBitVector :: (Num t, Eq t) =>
+      Vector t  -- ^Input bit-vector
+   -> Bool  -- ^Output boolean value
+isBitVector NullV = True
+isBitVector (x:>xs) = (x `elem` [0, 1]) && isBitVector xs
+
+-- |To transform the input integer to a bit-vector with specified number of
+-- bits.
+intToBitVector :: Int  -- ^Num of the bits
+       -> Integer  -- ^The input integer
+       -> BitVector -- ^The output bit-vector
+intToBitVector bits n | n >= 0 && n < 2^(bits-1) 
+       = intToBitVector' bits n
+intToBitVector bits n | n < 0 && abs n <= 2^(bits-1) 
+       = intToBitVector' bits (n + 2^bits)
+intToBitVector _    _ | otherwise = 
+      error "intToBitvector : Number out of range!" 
+
+-- |Helper function of 'intToBitVector'.
+intToBitVector' :: (Num a, Ord a1, Num a1, Integral t) => 
+       t -> a1 -> Vector a
+intToBitVector' 0    _ = NullV
+intToBitVector' bits n = if n >= 2^(bits-1) then
+         1 :> intToBitVector' (bits-1) (n - 2^(bits-1))
+       else  
+         0 :> intToBitVector' (bits-1) n
+
+-- |To transform the input bit-vecotr to an integer.
+bitVectorToInt :: BitVector -> Integer
+bitVectorToInt (1:>xv) | isBitVector xv 
+       = bitVectorToInt' xv (lengthV xv) - 2 ^ lengthV xv 
+bitVectorToInt (0:>xv) | isBitVector xv 
+       = bitVectorToInt' xv (lengthV xv)
+bitVectorToInt _ = error "bitVectorToInt: Vector is not a BitVector!"
+
+
+-- |Helper function of 'bitVectorToInt'.
+bitVectorToInt' :: (Integral a, Num t) => Vector t -> a -> t
+bitVectorToInt' NullV   _   = 0
+bitVectorToInt' (x:>xv) bit = x * 2^(bit-1) + bitVectorToInt' xv (bit-1)
+
+data Parity = Even | Odd deriving (Show, Eq)
+
+-- |To add even parity bit on the bit-vector in the tail.
+addEvenParityBit :: (Num a, Eq a) => Vector a -> Vector a
+addEvenParityBit = addParityBit Even
+-- |To add odd parity bit on the bit-vector in the tail.
+addOddParityBit :: (Num a, Eq a) => Vector a -> Vector a
+addOddParityBit  = addParityBit Odd
+
+addParityBit :: (Num a, Eq a) => Parity -> Vector a -> Vector a
+addParityBit p v 
+  | isBitVector v = case p of
+    Even -> resZero even
+    Odd -> resZero (not even)
+  | otherwise =  error "addParity: Vector is not a BitVector"  
+ where even = evenNumber v 
+       resZero b = v <+> unitV (if b then 0 else 1)
+
+
+-- |To remove the parity bit in the tail.
+removeParityBit :: (Num t, Eq t) => Vector t -> Vector t
+removeParityBit v 
+ | isBitVector v = takeV (lengthV v - 1) v
+ | otherwise = error "removeParityBit: Vector is not a BitVector "
+
+-- |To check the even parity of the bit-vector.
+isEvenParity :: (Num t, Eq t) => Vector t -> Bool
+isEvenParity = isParityCorrect Even
+
+-- |To check the odd parity of the bit-vector.
+isOddParity :: (Num t, Eq t) => Vector t -> Bool
+isOddParity = isParityCorrect Odd
+
+isParityCorrect :: (Num t, Eq t) => Parity -> Vector t -> Bool 
+isParityCorrect Even xv = evenNumber xv
+isParityCorrect Odd xv  = not $ evenNumber xv 
+
+evenNumber :: (Num t, Eq t) => Vector t -> Bool
+evenNumber NullV   = True
+evenNumber (0:>xv) = xor False (evenNumber xv)
+evenNumber (1:>xv) = xor True (evenNumber xv)
+evenNumber (_:>_) = error "evenNumber: Vector is not a BitVector "
+         
+xor :: Bool -> Bool -> Bool
+xor True  False = True
+xor False True  = True
+xor _     _     = False 
+
diff --git a/src/ForSyDe/Shallow/Utility/Matrix.hs b/src/ForSyDe/Shallow/Utility/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/ForSyDe/Shallow/Utility/Matrix.hs
@@ -0,0 +1,352 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  ForSyDe.Shallow.Utility.Matrix
+-- Copyright   :  (c) ForSyDe Group, KTH 2019
+-- License     :  BSD-style (see the file LICENSE)
+-- 
+-- Maintainer  :  forsyde-dev@kth.se
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- This module defines the data type 'Matrix' and corresponding
+-- functions. It is a shallow interpretation of 2D arrays and is used
+-- for quick prototyping of array algorithms and skeletons, although
+-- itself is a type synonym for a `Vector` of `Vector`s. Therefore
+-- this module is simply a collection of utility functions on 2D
+-- vectors used mainly for convenience. For a type-checked fixed-size
+-- data type for representing matrices, see the
+-- <http://hackage.haskell.org/package/matrix matrix> or
+-- <http://hackage.haskell.org/package/repa REPA> packages.
+--
+-- __OBS:__ The lengths in the API documentation for function arguments
+-- are not type-safe, but rather suggestions for usage in designing
+-- matrix algorithms or skeletons.
+-----------------------------------------------------------------------------
+module ForSyDe.Shallow.Utility.Matrix (
+  Matrix, prettyMat,
+  -- * Queries
+  nullMat, sizeMat, wellFormedMat,
+  -- * Generators
+  matrix, fromMatrix, unitMat, indexMat,
+  -- * Functional skeletons
+  mapMat, zipWithMat, zipWith3Mat, reduceMat, dotVecMat, dotMatMat,
+  -- * Selectors
+  atMat, takeMat, dropMat, cropMat, groupMat, stencilMat,
+  -- * Permutators
+  zipMat, unzipMat,
+  rotateMat, reverseMat, transposeMat, replaceMat
+  ) where
+
+import ForSyDe.Shallow.Core.Vector
+import Data.List (intercalate)
+
+-- | 'Matrix' is simply a type synonym for vector of vectors. This
+-- means that /any/ function on 'Vector' works also on 'Matrix'.
+type Matrix a = Vector (Vector a)
+
+-- | Prints out to the terminal a matrix in a readable format, where
+-- all elements are right-aligned and separated by a custom separator.
+--
+-- >>> let m = matrix 3 3 [1,2,3,3,100,4,12,32,67]
+-- >>> prettyMat "|" m
+--  1|  2| 3
+--  3|100| 4
+-- 12| 32|67
+prettyMat :: Show a
+          => String   -- ^ separator string
+          -> Matrix a -- ^ input matrix
+          -> IO ()
+prettyMat sep mat = mapM_ putStrLn $ fromVector $ printMat maxWdt strMat
+  where
+    maxWdt = reduceV (zipWithV max) $ mapMat length strMat
+    strMat = mapMat show mat
+    printMat w  = mapV (\row -> printRow w row)
+    printRow w  = intercalate sep . fromVector . zipWithV align w
+    align n str = replicate (n - length str) ' ' ++ str
+
+-- | Checks if a matrix is null. @<>@ and @<<>>@ are both null
+-- matrices.
+nullMat :: Matrix a -> Bool
+nullMat NullV = True
+nullMat (NullV:>NullV) = True
+nullMat _ = False
+
+-- | Returns the X and Y dimensions of matrix and checks if it is well formed.
+sizeMat :: Matrix a -> (Int,Int)
+sizeMat m = (x,y)
+  where
+    y = lengthV m
+    x = (lengthV . headV) (wellFormedMat m)
+
+-- | Checks if a matrix is well-formed, meaning that all its rows are
+-- of equal length. Returns the same matrix in case it is well-formed
+-- or throws an exception if it is ill-formed.
+wellFormedMat :: Matrix a -> Matrix a
+wellFormedMat NullV = error "matrix is null"
+wellFormedMat m@(_:>NullV) = m
+wellFormedMat m@(x:>xs)
+  | reduceV (&&) (mapV (\r -> lengthV r == lengthV x) xs) = m
+  | otherwise = error "matrix ill-formed: rows are of unequal lengths"
+
+groupEvery :: Int -> [a] -> [[a]]
+groupEvery _ [] = []
+groupEvery n l
+  | n < 0        = error $ "cannot group list by negative n: " ++ show n
+  | length l < n = error "input list cannot be split into all-equal parts"
+  | otherwise    = take n l : groupEvery n (drop n l)
+
+-- | Converts a list into a 'Matrix'. See example from 'prettyMat'.
+matrix :: Int      -- ^ number of columns (X dimension) @= x@
+       -> Int      -- ^ number of rows (Y dimension) @= y@
+       -> [a]      -- ^ list of values; /length/ = @x * y@
+       -> Matrix a -- ^ 'Matrix' of values; /size/ = @(x,y)@
+matrix x y = vector . map vector . groupEvery x . check
+  where
+    check l | length l == x * y = l
+            | otherwise
+      = error $ "cannot form matrix (" ++ show x ++ ","
+              ++ show y ++ ") from a list with "
+              ++ show (length l) ++ " elements"
+
+-- | Converts a matrix back to a list.
+fromMatrix :: Matrix a -- ^ /size/ = @(x,y)@
+           -> [a]      -- ^ /length/ = @x * y@
+fromMatrix = concatMap fromVector . fromVector
+
+-- | Creates a unit (i.e. singleton) matrix, which is a matrix with
+-- only one element.
+unitMat :: a -> Matrix a -- ^ /size/ = @(1,1)@
+unitMat a = (a:>NullV):>NullV
+
+-- | Returns an /infinite matrix/ with (X,Y) index pairs. You need to
+-- zip it against another (finite) matrix or to extract a finite
+-- subset in order to be useful (see example below).
+--
+-- >>> prettyMat " " $ takeMat 3 4 indexMat 
+-- (0,0) (1,0) (2,0)
+-- (0,1) (1,1) (2,1)
+-- (0,2) (1,2) (2,2)
+-- (0,3) (1,3) (2,3)
+indexMat :: Matrix (Int, Int)
+indexMat = zipWithMat (,) colix rowix
+  where
+    colix = vector $ repeat $ vector [0..]
+    rowix = transposeMat colix 
+
+-- | Maps a function on every value of a matrix.
+--
+-- __OBS:__ this function does not check if the output matrix is well-formed.
+mapMat :: (a -> b)
+       -> Matrix a -- ^ /size/ = @(xa,ya)@
+       -> Matrix b -- ^ /size/ = @(xa,ya)@
+mapMat = mapV . mapV
+
+-- | Applies a binary function pair-wise on each element in two matrices.
+--
+-- __OBS:__ this function does not check if the output matrix is well-formed.
+zipWithMat :: (a -> b -> c)
+           -> Matrix a -- ^ /size/ = @(xa,ya)@
+           -> Matrix b -- ^ /size/ = @(xb,yb)@
+           -> Matrix c -- ^ /size/ = @(minimum [xa,xb], minimum [ya,yb])@
+zipWithMat f = zipWithV (zipWithV f)
+
+-- | Applies a function 3-tuple-wise on each element in three matrices.
+--
+-- __OBS:__ this function does not check if the output matrix is well-formed.
+zipWith3Mat :: (a -> b -> c -> d)
+            -> Matrix a -- ^ /size/ = @(xa,ya)@
+            -> Matrix b -- ^ /size/ = @(xb,yb)@
+            -> Matrix c -- ^ /size/ = @(xc,yc)@
+            -> Matrix d -- ^ /size/ = @(minimum [xa,xb,xc], minimum [ya,yb,yc])@
+zipWith3Mat f = zipWith3V (\a b c -> zipWith3V f a b c)
+
+-- | Reduces all the elements of a matrix to one element based on a
+-- binary function.
+--
+-- >>> let m = matrix 3 3 [1,2,3,11,12,13,21,22,23]
+-- >>> reduceMat (+) m
+-- 108
+reduceMat :: (a -> a -> a) -> Matrix a -> a
+reduceMat f = reduceV f . mapV (reduceV f)
+
+-- | Pattern implementing the template for a dot operation between a
+-- vector and a matrix.
+--
+-- >>> let mA = matrix 4 4 [1,-1,1,1,  1,-1,-1,-1,  1,1,-1,1,  1,1,1,-1]
+-- >>> let y  = vector[1,0,0,0]
+-- >>> dotVecMat (+) (*) mA y
+-- <1,1,1,1>
+dotVecMat :: (a -> a -> a) -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
+          -> (b -> a -> a) -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
+          -> Matrix b      -- ^ /size/ = @(xa,ya)@
+          -> Vector a      -- ^ /length/ = @xa@
+          -> Vector a      -- ^ /length/ = @xa@
+dotVecMat f g mA y = mapV (\x -> reduceV f $ zipWithV g x y) mA
+
+-- | Pattern implementing the template for a dot operation between two
+-- matrices.
+--
+-- >>> let mA = matrix 4 4 [1,-1,1,1,  1,-1,-1,-1,  1,1,-1,1,  1,1,1,-1]
+-- >>> prettyMat " " $ dotMatMat (+) (*) mA mA
+-- 2 -2  2  2
+-- 2 -2 -2 -2
+-- 2  2  2 -2
+-- 2  2 -2  2
+dotMatMat :: (a -> a -> a) -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
+          -> (b -> a -> a) -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
+          -> Matrix b      -- ^ /size/ = @(xa,ya)@
+          -> Matrix a      -- ^ /size/ = @(ya,xa)@
+          -> Matrix a      -- ^ /size/ = @(xa,xa)@
+dotMatMat f g m = mapV (dotVecMat f g m) . transposeMat
+
+-- | Returns the element of a matrix at a certain position.
+--
+-- >>> let m = matrix 3 3 [1,2,3,11,12,13,21,22,23]
+-- >>> atMat 2 1 m
+-- 13
+atMat :: Int       -- ^ X index starting from zero
+      -> Int       -- ^ Y index starting from zero
+      -> Matrix a
+      -> a
+atMat x y mat = (mat `atV` y) `atV` x
+
+-- | Returns the upper-left part of a matrix until a specific
+-- position.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " $ takeMat 2 2 m
+--  1  2
+-- 11 12
+takeMat :: Int       -- ^ X index starting from zero
+        -> Int       -- ^ Y index starting from zero
+        -> Matrix a
+        -> Matrix a
+takeMat x y = mapV (takeV x) . takeV y
+
+-- | Returns the upper-left part of a matrix until a specific
+-- position.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " $ dropMat 2 2 m
+-- 23 24
+-- 33 34
+dropMat :: Int       -- ^ X index starting from zero
+        -> Int       -- ^ Y index starting from zero
+        -> Matrix a
+        -> Matrix a
+dropMat x y = mapV (dropV x) . dropV y
+
+-- | Crops a section of a matrix.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " m
+--  1  2  3  4
+-- 11 12 13 14
+-- 21 22 23 24
+-- 31 32 33 34
+-- >>> prettyMat " " $ cropMat 2 3 1 1 m
+-- 12 13
+-- 22 23
+-- 32 33
+cropMat :: Int      -- ^ crop width  = @w@
+        -> Int      -- ^ crop height = @h@
+        -> Int      -- ^ X start position = @x0@
+        -> Int      -- ^ Y start position = @y0@
+        -> Matrix a -- ^ /size/ = @(xa,ya)@
+        -> Matrix a -- ^ /size/ = @(minimum [w,xa-x0], minimum [h,xa-x0])@
+cropMat w h pX pY = takeMat w h . dropMat pX pY
+
+-- cropMat w h pX pY = mapV (crop w pX) . crop h pY
+--   where crop size pos = dropV pos . takeV (pos + size) 
+
+-- | Groups a matrix into smaller equallly-shaped matrices.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " $ groupMat 2 2 m
+--   <<1,2>,<11,12>>   <<3,4>,<13,14>>
+-- <<21,22>,<31,32>> <<23,24>,<33,34>>
+groupMat :: Int      -- ^ width of groups = @w@
+         -> Int      -- ^ height of groups = @h@
+         -> Matrix a -- ^ /size/ = @(xa,ya)@
+         -> Matrix (Matrix a) -- ^ /size/ = @(xa `div` w,ya `div` h)@
+groupMat w h = mapV transposeMat . groupV h . mapV (groupV w)
+
+
+-- | Returns a stencil of neighboring elements for each possible
+-- element in a vector.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " $ stencilMat 2 2 m
+--   <<1,2>,<11,12>>   <<2,3>,<12,13>>   <<3,4>,<13,14>>
+-- <<11,12>,<21,22>> <<12,13>,<22,23>> <<13,14>,<23,24>>
+-- <<21,22>,<31,32>> <<22,23>,<32,33>> <<23,24>,<33,34>>
+stencilMat :: Int -> Int -> Matrix a -> Matrix (Matrix a)
+stencilMat r c = arrange . groupCols . groupRows
+  where
+    groupRows =         mapV (takeV r) . dropFromEnd r . tailsV
+    groupCols = mapMat (mapV (takeV c) . dropFromEnd c . tailsV)
+    arrange   = mapV transposeMat
+    dropFromEnd n v = takeV (lengthV v - n) v
+
+-- | Reverses the order of elements in a matrix
+reverseMat :: Matrix a -> Matrix a
+reverseMat = reverseV . mapV reverseV
+
+-- | Pattern which "rotates" a matrix. The rotation is controled with
+-- the /x/ and /y/ index arguments as following:
+--
+-- * @(> 0)@ : rotates the matrix right/down with the corresponding
+-- number of positions.
+-- 
+-- * @(= 0)@ : does not modify the position for that axis.
+-- 
+-- * @(< 0)@ : rotates the matrix left/up with the corresponding
+-- number of positions.
+--
+-- >>> let m = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> prettyMat " " $ rotateMat (-1) 1 m
+-- 32 33 34 31
+--  2  3  4  1
+-- 12 13 14 11
+-- 22 23 24 21
+rotateMat :: Int -- ^ index on X axis
+          -> Int -- ^ index on Y axis
+          -> Vector (Vector a)
+          -> Vector (Vector a)
+rotateMat x y = rotateV y . mapV (rotateV x)
+
+-- | Transposes a matrx.
+transposeMat :: Matrix a -- ^ /size/ = @(x,y)@
+             -> Matrix a -- ^ /size/ = @(y,x)@
+transposeMat NullV = NullV
+transposeMat (NullV:>xss) = transposeMat xss
+transposeMat rows = (mapV headV rows) :> transposeMat (mapV tailV rows)
+
+-- | Replaces a part of matrix with another (smaller) part, starting
+-- from an arbitrary position.
+--
+-- >>> let m  = matrix 4 4 [1,2,3,4,11,12,13,14,21,22,23,24,31,32,33,34]
+-- >>> let m1 = matrix 2 2 [101,202,303,404]
+-- >>> prettyMat " " $ replaceMat 1 1 m1 m
+--  1   2   3  4
+-- 11 101 202 14
+-- 21 303 404 24
+-- 31  32  33 34
+replaceMat :: Int -> Int -> Matrix a -> Matrix a -> Matrix a
+replaceMat x y mask = replace y h (zipWithV (\m o -> replace x w (\_ -> m) o) mask)
+  where
+    (w,h) = sizeMat mask
+    replace start size replaceF vec
+      = let begin  = takeV start vec
+            middle = replaceF $ dropV start $ takeV (start + size) vec
+            end    = dropV (start + size) vec
+        in begin <+> middle <+> end
+
+zipMat :: Matrix a     -- ^ /size/ = @(xa,ya)@
+       -> Matrix b     -- ^ /size/ = @(xb,yb)@
+       -> Matrix (a,b) -- ^ /size/ = @(minimum [xa,xb], minimum [ya,yb])@
+zipMat = zipWithMat (,)
+
+unzipMat :: Matrix (a,b)         -- ^ /size/ = @(x,y)@
+         -> (Matrix a, Matrix b) -- ^ /size/ = @(x,y)@ and @(x,y)@
+unzipMat = unzipV . mapV unzipV
diff --git a/test/Props.hs b/test/Props.hs
new file mode 100644
--- /dev/null
+++ b/test/Props.hs
@@ -0,0 +1,71 @@
+import ForSyDe.Shallow
+import Test.QuickCheck
+
+type Rate = Positive Int
+
+instance (Arbitrary a) => Arbitrary (Signal a) where
+  arbitrary = do
+    x <- arbitrary
+    return (signal x)
+
+instance Arbitrary a => Arbitrary (AbstExt a) where
+  arbitrary = do
+    x <- arbitrary
+    return (Prst x)    
+
+countEvent _ NullS = 0
+countEvent a (x :- xs) | a == x    = 1 + countEvent a xs
+                       | otherwise = countEvent a xs   
+
+-- SY Process Properties
+
+prop_delaySY xs = 1 + lengthS xs == lengthS (delaySY 1 xs)
+  where types = xs :: Signal Int
+
+prop_mealySY xs = lengthS xs == lengthS (mealySY (+) (-) 1 xs)
+  where types = xs :: Signal Int
+
+prop_mapSY xs = lengthS xs == lengthS (mapSY (+1) xs)
+  where types = xs :: Signal Int
+
+prop_zipWithSY xs ys = min (lengthS xs) (lengthS ys) == lengthS (zipWithSY (+) xs ys)
+  where types = (xs :: Signal Int, ys :: Signal Int)
+
+prop_AbstSY1 xs = inputAbst xs == outputAbst xs
+  where types = xs :: Signal (AbstExt Int)
+        inputAbst    = countEvent Abst
+        outputAbst x = countEvent Abst (mapSY (psi (+1)) x)
+     
+
+-- SDF Process Properties
+
+prop_feedbackSDF xs = (lengthS xs) `div` (3 * 2) >= (lengthS out) `div` (2 * 3)
+  where types = xs :: Signal Int
+        out   = actor21SDF (3,1) 3 (\x y -> (+(head y)) <$> x) xs st
+        st    = delaySDF [1,1,1] $ actor21SDF (2,4) 2 (zipWith (+)) st xs
+
+prop_actor11SDF gc gp xs = rateI >= rateO && rateI <= rateO + 1
+  where types  = (gc :: Rate, gp :: Rate, xs :: Signal Int)
+        (c, p) = (getPositive gc, getPositive gc)
+        rateI  = lengthS xs  `div` c
+        rateO  = lengthS out `div` p
+        out    = actor11SDF c p (take p . repeat . head) xs
+
+prop_actor21SDF gc gp xs ys = rateI1 >= rateO && rateI2 >= rateO && (min rateI1 rateI2) <= rateO + 1
+  where types  = (gc :: (Rate,Rate), gp :: Rate, xs :: Signal Int, ys :: Signal Int)
+        (c1, c2, p) = (getPositive $ fst gc, getPositive $ snd gc, getPositive gp)
+        rateI1 = lengthS xs  `div` c1
+        rateI2 = lengthS ys  `div` c2
+        rateO  = lengthS out `div` p
+        out    = actor21SDF (c1,c2) p (\x -> take p . repeat . head) xs ys
+
+main = do
+  let runTest s prop = putStr (s ++ " ") >> quickCheck prop
+  runTest "SY delay num events" prop_delaySY
+  runTest "SY mealy num events" prop_mealySY
+  runTest "SY map num events" prop_mapSY
+  runTest "SY map num absents" prop_AbstSY1
+  runTest "SY zipWith num events" prop_zipWithSY
+  runTest "SDF feedback tokens" prop_zipWithSY
+  runTest "SDF actor11 tokens" prop_actor11SDF
+  runTest "SDF actor21 tokens" prop_actor21SDF
diff --git a/test/Spec.hs b/test/Spec.hs
deleted file mode 100644
--- a/test/Spec.hs
+++ /dev/null
@@ -1,62 +0,0 @@
-import ForSyDe.Shallow
-import Test.Hspec
-
--- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
-test_mulacc :: Integer -> Signal Integer
-test_mulacc n = sim_run $ mulacc constant1 siggen1
-  where
-    sim_run val = zipWithSY (\v _ -> v) val ticks 
-    mulacc a b  = result
-      where addi1  = comb2SY (*) a b         -- Multiplier
-            acci   = comb2SY (+) addi1 addi2 -- Adder
-            addi2  = delaySY 0 acci          -- Accumulator
-            result = acci                    -- Output of the system
-    constant1 = sourceSY id   3 :: Signal Integer
-    siggen1   = sourceSY (+1) 1 :: Signal Integer
-    ticks     = signal [1..n]   :: Signal Integer
-
-
--- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
-test_fibonacciRabbitsDeath :: Integer -> Signal Integer
-test_fibonacciRabbitsDeath months = fibonacciRabbitsDeath $ signal [1..months]
-  where
-    fibonacciRabbitsDeath ticks = zipWith4SY fusion n a d ticks
-      where n = newborns a
-            a = adults n
-            d = dead n
-            fusion x y z ctrl = x + y - z
-    newborns = delaySY 1
-    adults   = mooreSY nsf out 0
-      where nsf state input = (state + input)
-            out state = state
-    dead = delaynSY 0 4
-
--- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
-test_adaptiveAmp :: Integer -> Signal Integer
-test_adaptiveAmp n = sout
-    where s1   = p1 s3 sin -- Process p1
-          sout = p2 s1     -- Process p2
-          s2   = p3 sout   -- Process p3
-          s3   = p4 s2     -- Process p4
-          sin  = signal [10..n]
-          p1 = zipUs 1 5
-          p2 = mapU 1 mult
-            where mult [([control], signal)] = map (* control) signal
-          p3 = scanU (\_ -> 5) g 10  
-            where g :: (Ord a, Num a) => a -> [a] -> a
-                  g state signal 
-                    | sum signal > 500  = state - 1
-                    | sum signal < 400  = state + 1
-                    | otherwise         = state
-          p4 = initU [10] 
-
-main :: IO ()
-main = hspec $ do describe "ForSyDe.Shallow : " $ lab2tests
-  where
-    lab2tests = do
-      it "SY Multiply Accumulator" $ test_mulacc 10
-        `shouldBe`(read "{3,9,18,30,45,63,84,108,135,165}" :: Signal Integer)
-      it "SY Fibonacci Reproduction" $ test_fibonacciRabbitsDeath 10
-        `shouldBe`(read "{1,1,2,3,4,8,12,20,32,52}" :: Signal Integer)
-      it "U Adaptive Amplifier" $ test_adaptiveAmp 20
-        `shouldBe`(read "{100,110,120,130,140,135,144,153,162,171}" :: Signal Integer)
diff --git a/test/TestDocs.hs b/test/TestDocs.hs
new file mode 100644
--- /dev/null
+++ b/test/TestDocs.hs
@@ -0,0 +1,11 @@
+import Test.DocTest
+
+main = doctest
+  [ "-isrc"
+  , "src/ForSyDe/Shallow/MoC/Synchronous/Lib.hs"
+  , "src/ForSyDe/Shallow/MoC/SDF.hs"
+  , "src/ForSyDe/Shallow/MoC/CSDF.hs"
+  , "src/ForSyDe/Shallow/MoC/SADF.hs"
+  , "src/ForSyDe/Shallow/Core/Vector.hs"
+  , "src/ForSyDe/Shallow/Utility/Matrix.hs"
+  ]
diff --git a/test/Unit.hs b/test/Unit.hs
new file mode 100644
--- /dev/null
+++ b/test/Unit.hs
@@ -0,0 +1,126 @@
+import ForSyDe.Shallow
+import Test.Hspec
+
+-- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
+test_mulacc :: Integer -> Signal Integer
+test_mulacc n = sim_run $ mulacc constant1 siggen1
+  where
+    sim_run val = zipWithSY (\v _ -> v) val ticks 
+    mulacc a b  = result
+      where addi1  = comb2SY (*) a b         -- Multiplier
+            acci   = comb2SY (+) addi1 addi2 -- Adder
+            addi2  = delaySY 0 acci          -- Accumulator
+            result = acci                    -- Output of the system
+    constant1 = sourceSY id   3 :: Signal Integer
+    siggen1   = sourceSY (+1) 1 :: Signal Integer
+    ticks     = signal [1..n]   :: Signal Integer
+
+
+-- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
+test_fibonacciRabbitsDeath :: Integer -> Signal Integer
+test_fibonacciRabbitsDeath months = fibonacciRabbitsDeath $ signal [1..months]
+  where
+    fibonacciRabbitsDeath ticks = zipWith4SY fusion n a d ticks
+      where n = newborns a
+            a = adults n
+            d = dead n
+            fusion x y z ctrl = x + y - z
+    newborns = delaySY 1
+    adults   = mooreSY nsf out 0
+      where nsf state input = (state + input)
+            out state = state
+    dead = delaynSY 0 4
+
+-- | Taken from <https://github.com/forsyde/forsyde-shallow-examples>
+test_adaptiveAmp :: Integer -> Signal Integer
+test_adaptiveAmp n = sout
+    where s1   = p1 s3 sin -- Process p1
+          sout = p2 s1     -- Process p2
+          s2   = p3 sout   -- Process p3
+          s3   = p4 s2     -- Process p4
+          sin  = signal [10..n]
+          p1 = zipUs 1 5
+          p2 = mapU 1 mult
+            where mult [([control], signal)] = map (* control) signal
+          p3 = scanU (\_ -> 5) g 10  
+            where g :: (Ord a, Num a) => a -> [a] -> a
+                  g state signal 
+                    | sum signal > 500  = state - 1
+                    | sum signal < 400  = state + 1
+                    | otherwise         = state
+          p4 = initU [10] 
+
+-- | Testing feedback loops in in SDF library
+test_toySDF :: Signal Double -> Signal Double
+test_toySDF s_in = s_out where
+    s_1   = p_add s_in s_2
+    s_2   = p_delay s_1
+    s_out = p_average s_2   
+  -- Process specification
+    p_add s1 s2 = actor21SDF (1,1) 1 add s1 s2
+    p_delay s   = delaySDF [0] s
+    p_average s = actor11SDF 3 1 average s
+  -- Function definition
+    add [x] [y] = [x + y]
+    average [x1,x2,x3] = [(x1 + x2 + x3) / 3.0]
+s_test = signal [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
+
+-- test1: CSDF graph from the paper Cyclo-Static Dataflow
+test_toyCSDF :: Num a => Signal a
+test_toyCSDF = s3
+  where s3 = delayCSDF [1,1] s2
+        s2 = v2 s1
+        s1 = v1 s4
+        s4 = v3 s3
+        v1 = actor11CSDF [(1, 1, \[a] -> [a]), (1, 0, \_ -> []), (1, 0, \_ -> [])]
+        v2 = actor11CSDF [(1, 0, \_ -> []), (1, 2, \[a] -> [a, 2*a])]
+        v3 = actor11CSDF [(1, 3, \[a] -> [a, 2*a, 3*a])]
+
+test_antiWindUpSADF :: (Num a, Ord a) => Signal a -> Signal a
+test_antiWindUpSADF input = output
+  where
+    output = integrator c1 s1 s3
+    s3 = delaySADF [0] output
+    s1 = kernel11SADF c2 input
+    (c1, c2) = detector s3 input
+    integrator = kernel21SADF
+    -- Detector
+    detector :: (Num a, Ord a) => Signal a -> Signal a
+             -> (Signal ((Int, Int), Int, [a] -> [a] -> [a]),
+                 Signal (Int, Int, [a] -> [a]))
+    detector = detector22SADF (1,1) f g 1
+    -- State transition function for the detector
+    f :: (Num a, Ord a) => Int -> [a] -> [a] -> Int
+    f 1 [y] [v] = if (y > 100 && v > 0 || y < (-100) && v < 0)
+                  then 2
+                  else 1
+    f 2 [y] [v] = if (y > 100 && v > 0 || y < (-100) && v < 0)
+                  then 2
+                  else 1
+    -- Output function for the detector
+    g :: Num a => Int -> ((Int, Int), ([((Int, Int), Int, [a] -> [a] -> [a])], [(Int, Int, [a] -> [a])]))
+    g 1 = ((1,1),
+           ([((1,1), 1, \[a] [b] -> [a+b])],
+            [(1, 1, \[a] -> [a])]))
+    g 2 = ((1,1),
+           ([((0,1), 1, \_ [b] -> [b])],
+            [(1, 0, \[a] -> [])]))
+s2_test = signal $ [10..110]
+
+                 
+main :: IO ()
+main = hspec $ do describe "ForSyDe.Shallow : " $ lab2tests
+  where
+    lab2tests = do
+      it "SY Multiply Accumulator" $ test_mulacc 10
+        `shouldBe`(read "{3,9,18,30,45,63,84,108,135,165}" :: Signal Integer)
+      it "SY Fibonacci Reproduction" $ test_fibonacciRabbitsDeath 10
+        `shouldBe`(read "{1,1,2,3,4,8,12,20,32,52}" :: Signal Integer)
+      it "U Adaptive Amplifier" $ test_adaptiveAmp 20
+        `shouldBe`(read "{100,110,120,130,140,135,144,153,162,171}" :: Signal Integer)
+      it "SDF Feedback Loop" $ test_toySDF s_test
+        `shouldBe`(read "{1.3333333333333333,10.333333333333334}" :: Signal Double)
+      it "CSDF Feedback Loop" $ takeS 10 test_toyCSDF
+        `shouldBe`(read "{1,1,1,2,2,4,4,8,8,16}" :: Signal Integer)
+      it "SADF Anti Wind-up System" $ takeS 10 (test_antiWindUpSADF s2_test)
+        `shouldBe`(read "{10,21,33,46,60,75,91,108,108,108}" :: Signal Integer)
