diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 2011 Jakub Waszczuk, 2012 IPI PAN
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/crf-chain2-tiers.cabal b/crf-chain2-tiers.cabal
new file mode 100644
--- /dev/null
+++ b/crf-chain2-tiers.cabal
@@ -0,0 +1,57 @@
+name:               crf-chain2-tiers
+version:            0.1.0
+synopsis:           Second-order, tiered, constrained, linear conditional random fields
+description:
+    The library provides implementation of the second-order, linear
+    conditional random fields (CRFs) with position-wise constraints
+    imposed over label values.  Each label consists of a vector of
+    smaller, atomic labels, and over each tier (layer) a separate
+    set of model features is defined.
+license:            BSD3
+license-file:       LICENSE
+cabal-version:      >= 1.6
+copyright:          Copyright (c) 2013 IPI PAN
+author:             Jakub Waszczuk
+maintainer:         waszczuk.kuba@gmail.com
+stability:          experimental
+category:           Natural Language Processing
+homepage:           https://github.com/kawu/crf-chain2-tiers
+build-type:         Simple
+
+library
+    hs-source-dirs: src
+
+    build-depends:
+        base >= 4 && < 5
+      , containers
+      , array
+      , vector
+      , binary
+      , vector-binary
+      , monad-codec >= 0.2 && < 0.3
+      , data-lens
+      , comonad-transformers
+      , logfloat
+      , parallel
+      , sgd >= 0.3 && < 0.4
+
+    exposed-modules:
+        Data.CRF.Chain2.Tiers
+      , Data.CRF.Chain2.Tiers.Dataset.Internal
+      , Data.CRF.Chain2.Tiers.Dataset.External
+      , Data.CRF.Chain2.Tiers.Dataset.Codec
+      , Data.CRF.Chain2.Tiers.Feature
+      , Data.CRF.Chain2.Tiers.Model
+      , Data.CRF.Chain2.Tiers.Inference
+      , Data.CRF.Chain2.Tiers.Array
+
+    other-modules:
+        Data.CRF.Chain2.Tiers.Util
+      , Data.CRF.Chain2.Tiers.DP
+        
+
+    ghc-options: -Wall -O2
+
+source-repository head
+    type: git
+    location: https://github.com/kawu/crf-chain2-tiers.git
diff --git a/src/Data/CRF/Chain2/Tiers.hs b/src/Data/CRF/Chain2/Tiers.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers.hs
@@ -0,0 +1,186 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE PatternGuards #-}
+
+
+module Data.CRF.Chain2.Tiers
+( 
+-- * CRF
+  CRF (..)
+, size
+, prune
+
+-- * Training
+, train
+, reTrain
+
+-- * Tagging
+, tag
+
+-- * Modules
+, module Data.CRF.Chain2.Tiers.Dataset.External
+, module Data.CRF.Chain2.Tiers.Feature
+) where
+
+
+import           System.IO (hSetBuffering, stdout, BufferMode (..))
+import           Control.Applicative ((<$>), (<*>))
+import           Data.Maybe (maybeToList)
+import qualified Data.Map as M
+import           Data.Binary (Binary, get, put)
+import qualified Data.Number.LogFloat as LogFloat
+import qualified Numeric.SGD as SGD
+import qualified Numeric.SGD.LogSigned as L
+
+import           Data.CRF.Chain2.Tiers.Dataset.Internal
+import           Data.CRF.Chain2.Tiers.Dataset.Codec
+import           Data.CRF.Chain2.Tiers.Dataset.External
+import           Data.CRF.Chain2.Tiers.Feature
+import           Data.CRF.Chain2.Tiers.Model
+import qualified Data.CRF.Chain2.Tiers.Inference as I
+
+
+----------------------------------------------------
+-- CRF model
+----------------------------------------------------
+
+
+-- | CRF model data.
+data CRF a b = CRF
+    { numOfLayers   :: Int
+    , codec         :: Codec a b
+    , model         :: Model }
+
+
+instance (Ord a, Ord b, Binary a, Binary b) => Binary (CRF a b) where
+    put CRF{..} = put numOfLayers >> put codec >> put model
+    get = CRF <$> get <*> get <*> get
+
+
+-- | Compute size (number of features) of the model.
+size :: CRF a b -> Int
+size CRF{..} = M.size (toMap model)
+
+
+-- | Discard model features with absolute values (in log-domain)
+-- lower than the given threshold.
+prune :: Double -> CRF a b -> CRF a b
+prune x crf =  crf { model = newModel } where
+    newModel = fromMap . M.fromList $
+        [ (feat, val)
+        | (feat, val) <- M.toList $ toMap (model crf)
+        , abs (LogFloat.logFromLogFloat val) > x ]
+
+
+----------------------------------------------------
+-- Training
+----------------------------------------------------
+
+
+-- | Train the CRF using the stochastic gradient descent method.
+train
+    :: (Ord a, Ord b)
+    => Int                          -- ^ Number of layers (tiers)
+    -> FeatSel                      -- ^ Feature selection
+    -> SGD.SgdArgs                  -- ^ SGD parameters
+    -> Bool                         -- ^ Store dataset on a disk
+    -> IO [SentL a b]               -- ^ Training data 'IO' action
+    -> IO [SentL a b]               -- ^ Evaluation data
+    -> IO (CRF a b)                 -- ^ Resulting model
+train numOfLayers featSel sgdArgs onDisk trainIO evalIO = do
+    hSetBuffering stdout NoBuffering
+
+    -- Create codec and encode the training dataset
+    codec <- mkCodec numOfLayers    <$> trainIO
+    trainData_ <- encodeDataL codec <$> trainIO
+    SGD.withData onDisk trainData_ $ \trainData -> do
+
+    -- Encode the evaluation dataset
+    evalData_ <- encodeDataL codec <$> evalIO
+    SGD.withData onDisk evalData_ $ \evalData -> do
+
+    -- Train the model
+    model <- mkModel featSel <$> SGD.loadData trainData
+    para  <- SGD.sgd sgdArgs
+        (notify sgdArgs model trainData evalData)
+        (gradOn model) trainData (values model)
+    return $ CRF numOfLayers codec model { values = para }
+
+
+-- | Re-train the CRF using the stochastic gradient descent method.
+reTrain
+    :: (Ord a, Ord b)
+    => CRF a b                      -- ^ Existing CRF model
+    -> SGD.SgdArgs                  -- ^ SGD parameters
+    -> Bool                         -- ^ Store dataset on a disk
+    -> IO [SentL a b]               -- ^ Training data 'IO' action
+    -> IO [SentL a b]               -- ^ Evaluation data
+    -> IO (CRF a b)                 -- ^ Resulting model
+reTrain crf sgdArgs onDisk trainIO evalIO = do
+    hSetBuffering stdout NoBuffering
+
+    -- Encode the training dataset
+    trainData_ <- encodeDataL (codec crf) <$> trainIO
+    SGD.withData onDisk trainData_ $ \trainData -> do
+
+    -- Encode the evaluation dataset
+    evalData_ <- encodeDataL (codec crf) <$> evalIO
+    SGD.withData onDisk evalData_ $ \evalData -> do
+
+    -- Train the model
+    let model' = model crf
+    para  <- SGD.sgd sgdArgs
+        (notify sgdArgs model' trainData evalData)
+        (gradOn model') trainData (values model')
+    return $ crf { model = model' { values = para } }
+
+
+-- | Compute gradient on a dataset element.
+gradOn :: Model -> SGD.Para -> (Xs, Ys) -> SGD.Grad
+gradOn model para (xs, ys) = SGD.fromLogList $
+    [ (unFeatIx ix, L.fromPos val)
+    | (ft, val) <- presentFeats xs ys
+    , ix <- maybeToList (index curr ft) ] ++
+    [ (unFeatIx ix, L.fromNeg val)
+    | (ft, val) <- I.expectedFeatures curr xs
+    , ix <- maybeToList (index curr ft) ]
+  where
+    curr = model { values = para }
+
+
+notify
+    :: SGD.SgdArgs -> Model
+    -> SGD.Dataset (Xs, Ys)         -- ^ Training dataset
+    -> SGD.Dataset (Xs, Ys)         -- ^ Evaluaion dataset
+    -> SGD.Para -> Int -> IO ()
+notify SGD.SgdArgs{..} model trainData evalData para k
+    | doneTotal k == doneTotal (k - 1) = putStr "."
+    | SGD.size evalData > 0 = do
+        x <- I.accuracy (model { values = para }) <$> SGD.loadData evalData
+        putStrLn ("\n" ++ "[" ++ show (doneTotal k) ++ "] f = " ++ show x)
+    | otherwise =
+        putStrLn ("\n" ++ "[" ++ show (doneTotal k) ++ "] f = #")
+  where
+    doneTotal :: Int -> Int
+    doneTotal = floor . done
+    done :: Int -> Double
+    done i
+        = fromIntegral (i * batchSize)
+        / fromIntegral trainSize
+    trainSize = SGD.size trainData
+
+
+----------------------------------------------------
+-- Tagging
+----------------------------------------------------
+
+
+-- | Find the most probable label sequence.
+tag :: (Ord a, Ord b) => CRF a b -> Sent a b -> [[b]]
+tag CRF{..} sent
+    = onWords . decodeLabels codec
+    . I.tag model . encodeSent codec
+    $ sent
+  where
+    onWords xs =
+        [ unJust codec word x
+        | (word, x) <- zip sent xs ]
diff --git a/src/Data/CRF/Chain2/Tiers/Array.hs b/src/Data/CRF/Chain2/Tiers/Array.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Array.hs
@@ -0,0 +1,113 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE TupleSections #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+
+module Data.CRF.Chain2.Tiers.Array
+(
+-- * Array
+  Array
+, mkArray
+, unArray
+, (!?)
+
+-- * Bounds
+, Bounds (..)
+) where
+
+
+import           Control.Applicative ((<$>), (<*>))
+import           Control.Arrow (first)
+import           Data.Ix
+import           Data.Int (Int16)
+import           Data.Maybe (catMaybes)
+import           Data.List (foldl1')
+import qualified Data.Vector.Unboxed as U
+import           Data.Binary (Binary, get, put)
+import           Data.Vector.Binary ()
+
+
+--------------------------------
+-- Array
+--------------------------------
+
+
+-- | An unboxed array implemented in terms of an unboxed vector.  
+data Array i a = Array
+    { bounds    :: (i, i)
+    , array     :: U.Vector a }
+
+
+instance (Binary i, Binary a, U.Unbox a) => Binary (Array i a) where
+    put Array{..} = put bounds >> put array
+    get = Array <$> get <*> get
+
+
+-- | Construct array with a default dummy value.
+mkArray :: (Bounds i, U.Unbox a) => a -> [(i, a)] -> Array i a
+mkArray dummy xs = Array
+    { bounds    = (p, q)
+    , array     = zeroed U.// map (first ix) xs }
+  where
+    p       = foldl1' lower (map fst xs)
+    q       = foldl1' upper (map fst xs)
+    ix      = index (p, q)
+    size    = rangeSize (p, q)
+    zeroed  = U.replicate size dummy
+{-# INLINE mkArray #-}
+
+
+-- | Deconstruct the array.
+unArray :: (Bounds i, U.Unbox a) => Array i a -> [(i, a)]
+unArray ar = catMaybes
+    [ (i,) <$> (ar !? i)
+    | i <- range (bounds ar) ]
+{-# INLINE unArray #-}
+
+
+(!?) :: (Ix i, U.Unbox a) => Array i a -> i -> Maybe a
+Array{..} !? x = if inRange bounds x
+    -- TODO: Use unsafe indexing.
+    then Just (array U.! index bounds x)
+    else Nothing
+{-# INLINE (!?) #-}
+
+
+--------------------------------
+-- Bounds
+--------------------------------
+
+
+-- | An extended Ix class.
+class Ix i => Bounds i where
+    -- | A lower bound of two values.
+    lower :: i -> i -> i
+    -- | An upper bound of two values.
+    upper :: i -> i -> i
+
+
+instance Bounds Int16 where
+    lower x y = min x y
+    upper x y = max x y
+
+
+instance Bounds i => Bounds (i, i) where
+    lower (!x1, !y1) (!x2, !y2) =
+        ( lower x1 x2
+        , lower y1 y2 )
+    upper (!x1, !y1) (!x2, !y2) =
+        ( upper x1 x2
+        , upper y1 y2 )
+
+
+instance Bounds i => Bounds (i, i, i) where
+    lower (!x1, !y1, !z1) (!x2, !y2, !z2) =
+        ( lower x1 x2
+        , lower y1 y2
+        , lower z1 z2 )
+    upper (!x1, !y1, !z1) (!x2, !y2, !z2) =
+        ( upper x1 x2
+        , upper y1 y2
+        , upper z1 z2 )
diff --git a/src/Data/CRF/Chain2/Tiers/DP.hs b/src/Data/CRF/Chain2/Tiers/DP.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/DP.hs
@@ -0,0 +1,43 @@
+module Data.CRF.Chain2.Tiers.DP
+( table
+, flexible2
+, flexible3
+) where
+
+import qualified Data.Array as A
+import Data.Array ((!))
+import Data.Ix (range)
+
+table :: A.Ix i => (i, i) -> ((i -> e) -> i -> e) -> A.Array i e
+table bounds f = table' where
+    table' = A.listArray bounds
+           $ map (f (table' !)) 
+           $ range bounds
+
+down1 :: A.Ix i => (i, i) -> (i -> e) -> i -> e
+down1 bounds f = (!) down' where
+    down' = A.listArray bounds
+          $ map f
+          $ range bounds
+
+down2 :: (A.Ix i, A.Ix j) => (j, j) -> (j -> (i, i)) -> (j -> i -> e)
+      -> j -> i -> e
+down2 bounds1 bounds2 f = (!) down' where
+    down' = A.listArray bounds1
+        [ down1 (bounds2 i) (f i)
+        | i <- range bounds1 ]
+
+flexible2 :: (A.Ix i, A.Ix j) => (j, j) -> (j -> (i, i))  
+          -> ((j -> i -> e) -> j -> i -> e) -> j -> i -> e
+flexible2 bounds1 bounds2 f = (!) flex where
+    flex = A.listArray bounds1
+        [ down1 (bounds2 i) (f (flex !) i)
+        | i <- range bounds1 ]
+
+flexible3 :: (A.Ix j, A.Ix i, A.Ix k) => (k, k) -> (k -> (j, j))
+          -> (k -> j -> (i, i)) -> ((k -> j -> i -> e) -> k -> j -> i -> e)
+           -> k -> j -> i -> e
+flexible3 bounds1 bounds2 bounds3 f = (!) flex where
+    flex = A.listArray bounds1
+        [ down2 (bounds2 i) (bounds3 i) (f (flex !) i)
+        | i <- range bounds1 ]
diff --git a/src/Data/CRF/Chain2/Tiers/Dataset/Codec.hs b/src/Data/CRF/Chain2/Tiers/Dataset/Codec.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Dataset/Codec.hs
@@ -0,0 +1,258 @@
+module Data.CRF.Chain2.Tiers.Dataset.Codec
+( Codec
+, CodecM
+, obMax
+, lbMax
+
+, encodeWord'Cu
+, encodeWord'Cn
+, encodeSent'Cu
+, encodeSent'Cn
+, encodeSent
+
+, encodeWordL'Cu
+, encodeWordL'Cn
+, encodeSentL'Cu
+, encodeSentL'Cn
+, encodeSentL
+
+-- , encodeLabels
+, decodeLabel
+, decodeLabels
+
+, mkCodec
+, encodeData
+, encodeDataL
+, unJust
+) where
+
+
+import Control.Applicative ((<$>), (<*>), pure)
+import Control.Comonad.Trans.Store (store)
+import Data.Maybe (catMaybes, fromJust)
+import Data.Lens.Common (Lens(..), fstLens)
+import qualified Data.Set as S
+import qualified Data.Map as M
+import qualified Data.Vector as V
+import qualified Control.Monad.Codec as C
+
+import Data.CRF.Chain2.Tiers.Dataset.Internal
+import Data.CRF.Chain2.Tiers.Dataset.External
+
+
+-- | Codec internal data.  The first component is used to
+-- encode observations of type a, the second one is used to
+-- encode labels of type [b].
+type Codec a b =
+    ( C.AtomCodec a
+    , V.Vector (C.AtomCodec (Maybe b)) )
+
+
+-- | The maximum internal observation included in the codec.
+obMax :: Codec a b -> Ob
+obMax =
+    let idMax m = M.size m - 1
+    in  mkOb . idMax . C.to . fst
+
+
+-- | The maximum internal labels included in the codec.
+lbMax :: Codec a b -> [Lb]
+lbMax =
+    let idMax m = M.size m - 1
+    in  map (mkLb . idMax . C.to) . V.toList . snd
+
+
+obLens :: Lens (a, b) a
+obLens = fstLens
+
+
+lbLens :: Int -> Lens (a, V.Vector b) b
+lbLens k = Lens $ \(a, b) -> store
+    (\x -> (a, b V.// [(k, x)]))
+    (b V.! k)
+
+
+--------------------------------------
+-- Core
+--------------------------------------
+
+
+-- | The empty codec.  The label parts are initialized with Nothing
+-- members, which represent unknown labels.  It is taken in the model
+-- implementation on account because it is assigned to the
+-- lowest label code and the model assumes that the set of labels
+-- is of the {0, ..., 'lbMax'} form.
+--
+-- Codec depends on the number of layers.
+empty :: Ord b => Int -> Codec a b
+empty n =
+    let withNo = C.execCodec C.empty (C.encode C.idLens Nothing)
+    in  (C.empty, V.replicate n withNo)
+
+
+-- | Type synonym for the codec monad.
+type CodecM a b c = C.Codec (Codec a b) c
+
+
+-- | Encode the observation and update the codec (only in the encoding
+-- direction).
+encodeObU :: Ord a => a -> CodecM a b Ob
+encodeObU = fmap mkOb . C.encode' obLens
+
+
+-- | Encode the observation and do *not* update the codec.
+encodeObN :: Ord a => a -> CodecM a b (Maybe Ob)
+encodeObN = fmap (fmap mkOb) . C.maybeEncode obLens
+
+
+-- | Encode the label and update the codec.
+encodeLbU :: Ord b => [b] -> CodecM a b Cb
+encodeLbU xs = mkCb <$> sequence
+    [ mkLb <$> C.encode (lbLens k) (Just x)
+    | (x, k) <- zip xs [0..] ]
+
+
+-- | Encode the label and do *not* update the codec.
+-- In case the label is not a member of the codec,
+-- return the label code assigned to Nothing label.
+encodeLbN :: Ord b => [b] -> CodecM a b Cb
+encodeLbN xs =
+    let encode lens x = C.maybeEncode lens (Just x) >>= \mx -> case mx of
+            Just x' -> return x'
+            Nothing -> fromJust <$> C.maybeEncode lens Nothing
+    in  mkCb <$> sequence
+            [ mkLb <$> encode (lbLens k) x
+            | (x, k) <- zip xs [0..] ]
+
+
+-- | Decode the label within the codec monad.
+decodeLbC :: Ord b => Cb -> CodecM a b (Maybe [b])
+decodeLbC xs = sequence <$> sequence
+    [ C.decode (lbLens k) (unLb x)
+    | (x, k) <- zip (unCb xs) [0..] ]
+
+
+-- | Is label a member of the codec?
+hasLabel :: Ord b => Codec a b -> [b] -> Bool
+hasLabel cdc xs = and
+        [ M.member
+            (Just x)
+            (C.to $ snd cdc V.! k)
+        | (x, k) <- zip xs [0..] ]
+
+
+--------------------------------------
+-- On top of core
+--------------------------------------
+
+
+-- | Encode the labeled word and update the codec.
+encodeWordL'Cu :: (Ord a, Ord b) => WordL a b -> CodecM a b (X, Y)
+encodeWordL'Cu (word, choice) = do
+    x' <- mapM encodeObU (S.toList (obs word))
+    r' <- mapM encodeLbU (S.toList (lbs word))
+    let x = mkX x' r'
+    y  <- mkY <$> sequence
+    	[ (,) <$> encodeLbU lb <*> pure pr
+	| (lb, pr) <- (M.toList . unProb) choice ]
+    return (x, y)
+
+-- | Encodec the labeled word and do *not* update the codec.
+encodeWordL'Cn :: (Ord a, Ord b) => WordL a b -> CodecM a b (X, Y)
+encodeWordL'Cn (word, choice) = do
+    x' <- catMaybes <$> mapM encodeObN (S.toList (obs word))
+    r' <- mapM encodeLbN (S.toList (lbs word))
+    let x = mkX x' r'
+    y  <- mkY <$> sequence
+    	[ (,) <$> encodeLbN lb <*> pure pr
+	| (lb, pr) <- (M.toList . unProb) choice ]
+    return (x, y)
+
+-- | Encode the word and update the codec.
+encodeWord'Cu :: (Ord a, Ord b) => Word a b -> CodecM a b X
+encodeWord'Cu word = do
+    x' <- mapM encodeObU (S.toList (obs word))
+    r' <- mapM encodeLbU (S.toList (lbs word))
+    return $ mkX x' r'
+
+-- | Encode the word and do *not* update the codec.
+encodeWord'Cn :: (Ord a, Ord b) => Word a b -> CodecM a b X
+encodeWord'Cn word = do
+    x' <- catMaybes <$> mapM encodeObN (S.toList (obs word))
+    r' <- mapM encodeLbN (S.toList (lbs word))
+    return $ mkX x' r'
+
+-- | Encode the labeled sentence and update the codec.
+encodeSentL'Cu :: (Ord a, Ord b) => SentL a b -> CodecM a b (Xs, Ys)
+encodeSentL'Cu sent = do
+    ps <- mapM (encodeWordL'Cu) sent
+    return (V.fromList (map fst ps), V.fromList (map snd ps))
+
+-- | Encode the labeled sentence and do *not* update the codec.
+-- Substitute the default label for any label not present in the codec.
+encodeSentL'Cn :: (Ord a, Ord b) => SentL a b -> CodecM a b (Xs, Ys)
+encodeSentL'Cn sent = do
+    ps <- mapM (encodeWordL'Cn) sent
+    return (V.fromList (map fst ps), V.fromList (map snd ps))
+
+-- -- | Encode labels into an ascending vector of distinct label codes.
+-- encodeLabels :: Ord b => Codec a b -> [b] -> AVec Lb
+-- encodeLabels codec = mkAVec . C.evalCodec codec . mapM encodeLbN
+
+-- | Encode the labeled sentence with the given codec.  Substitute the
+-- default label for any label not present in the codec.
+encodeSentL :: (Ord a, Ord b) => Codec a b -> SentL a b -> (Xs, Ys)
+encodeSentL codec = C.evalCodec codec . encodeSentL'Cn
+
+-- | Encode the sentence and update the codec.
+encodeSent'Cu :: (Ord a, Ord b) => Sent a b -> CodecM a b Xs
+encodeSent'Cu = fmap V.fromList . mapM encodeWord'Cu
+
+-- | Encode the sentence and do *not* update the codec.
+encodeSent'Cn :: (Ord a, Ord b) => Sent a b -> CodecM a b Xs
+encodeSent'Cn = fmap V.fromList . mapM encodeWord'Cn
+
+-- | Encode the sentence using the given codec.
+encodeSent :: (Ord a, Ord b) => Codec a b -> Sent a b -> Xs
+encodeSent codec = C.evalCodec codec . encodeSent'Cn
+
+-- | Create codec on the basis of the labeled dataset.
+mkCodec :: (Ord a, Ord b) => Int -> [SentL a b] -> Codec a b
+mkCodec n = C.execCodec (empty n) . mapM_ encodeSentL'Cu
+
+-- | Encode the labeled dataset using the codec.  Substitute the default
+-- label for any label not present in the codec.
+encodeDataL :: (Ord a, Ord b) => Codec a b -> [SentL a b] -> [(Xs, Ys)]
+encodeDataL = map . encodeSentL
+
+-- | Encode the dataset with the codec.
+encodeData :: (Ord a, Ord b) => Codec a b -> [Sent a b] -> [Xs]
+encodeData = map . encodeSent
+
+-- | Decode the label.
+decodeLabel :: Ord b => Codec a b -> Cb -> Maybe [b]
+decodeLabel codec = C.evalCodec codec . decodeLbC
+
+-- | Decode the sequence of labels.
+decodeLabels :: Ord b => Codec a b -> [Cb] -> [Maybe [b]]
+decodeLabels codec = C.evalCodec codec . mapM decodeLbC
+
+-- | Return the label when 'Just' or one of the unknown values
+-- when 'Nothing'.
+unJust :: Ord b => Codec a b -> Word a b -> Maybe [b] -> [b]
+unJust _ _ (Just x) = x
+unJust codec word Nothing = case allUnk of
+    (x:_)   -> x
+    []      -> error "unJust: Nothing and all values known"
+  where
+    allUnk = filter (not . hasLabel codec) (S.toList $ lbs word)
+
+-- -- | Replace 'Nothing' labels with all unknown labels from
+-- -- the set of potential interpretations.
+-- unJusts :: Ord b => Codec a b -> Word a b -> [Maybe b] -> [b]
+-- unJusts codec word xs =
+--     concatMap deJust xs
+--   where
+--     allUnk = filter (not . hasLabel codec) (S.toList $ lbs word)
+--     deJust (Just x) = [x]
+--     deJust Nothing  = allUnk
diff --git a/src/Data/CRF/Chain2/Tiers/Dataset/External.hs b/src/Data/CRF/Chain2/Tiers/Dataset/External.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Dataset/External.hs
@@ -0,0 +1,59 @@
+-- | External data representation.
+
+module Data.CRF.Chain2.Tiers.Dataset.External
+( Word (obs, lbs)
+, mkWord
+, Sent
+, Prob (unProb)
+, mkProb
+, WordL
+, SentL
+) where
+
+import qualified Data.Set as S
+import qualified Data.Map as M
+
+-- | A word consists of a set of observations and a set of potential labels.
+data Word a b = Word {
+    -- | Set of observations.
+      obs   :: S.Set a
+    -- | Non-empty set of potential labels.
+    , lbs   :: S.Set [b] }
+    deriving (Show, Eq, Ord)
+
+-- | A word constructor which checks non-emptiness of the potential
+-- set of labels.
+mkWord :: S.Set a -> S.Set [b] -> Word a b
+mkWord _obs _lbs
+    | S.null _lbs   = error "mkWord: empty set of potential labels"
+    | otherwise     = Word _obs _lbs
+
+-- | A sentence of words.
+type Sent a b = [Word a b]
+
+-- | A probability distribution defined over elements of type a.
+-- All elements not included in the map have probability equal
+-- to 0.
+newtype Prob a = Prob { unProb :: M.Map a Double }
+    deriving (Show, Eq, Ord)
+
+-- | Construct the probability distribution.
+mkProb :: Ord a => [(a, Double)] -> Prob a
+mkProb =
+    Prob . normalize . M.fromListWith (+) . filter ((>0).snd)
+  where
+    normalize dist 
+        | M.null dist =
+            error "mkProb: no elements with positive probability"
+        | otherwise   =
+            let z = sum (M.elems dist)
+            in  fmap (/z) dist
+
+-- | A WordL is a labeled word, i.e. a word with probability distribution
+-- defined over labels.  We assume that every label from the distribution
+-- domain is a member of the set of potential labels corresponding to the
+-- word.  TODO: Ensure the assumption using the smart constructor.
+type WordL a b = (Word a b, Prob [b])
+
+-- | A sentence of labeled words.
+type SentL a b = [WordL a b]
diff --git a/src/Data/CRF/Chain2/Tiers/Dataset/Internal.hs b/src/Data/CRF/Chain2/Tiers/Dataset/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Dataset/Internal.hs
@@ -0,0 +1,279 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RecordWildCards #-}
+
+
+-- | Internal core data types.
+
+
+module Data.CRF.Chain2.Tiers.Dataset.Internal
+(
+-- * Basic types
+  Ob (..)
+, mkOb, unOb
+, Lb (..)
+, mkLb, unLb
+, FeatIx (..)
+, mkFeatIx, unFeatIx
+, CbIx
+
+-- * Complex label
+, Cb (..)
+, mkCb
+, unCb
+
+-- * Input element (word)
+, X (_unX, _unR)
+, Xs
+, mkX
+, unX
+, unR
+
+-- * Output element (choice)
+, Y (_unY)
+, Ys
+, mkY
+, unY
+
+-- * Indexing
+, lbAt
+, lbOn
+, lbNum
+, lbIxs
+) where
+
+
+import           Data.Binary (Binary, put, get)
+import           Data.Ix (Ix)
+import           Control.Applicative ((<$>), (<*>))
+import           Control.Arrow (second)
+import           Data.Int (Int16, Int32)
+import qualified Data.Array.Unboxed as A
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic.Base as G
+import qualified Data.Vector.Generic.Mutable as G
+import qualified Data.Number.LogFloat as L
+-- import qualified Data.Primitive.ByteArray as BA
+
+import           Data.CRF.Chain2.Tiers.Array (Bounds)
+
+----------------------------------------------------------------
+-- Basic types
+----------------------------------------------------------------
+
+
+-- | An observation.
+newtype Ob = Ob { _unOb :: Int32 }
+    deriving ( Show, Eq, Ord, Binary, A.IArray A.UArray
+             , G.Vector U.Vector, G.MVector U.MVector, U.Unbox )
+
+
+-- | Smart observation constructor.
+mkOb :: Int -> Ob
+mkOb = Ob . fromIntegral
+{-# INLINE mkOb #-}
+
+
+-- | Deconstract observation.
+unOb :: Ob -> Int
+unOb = fromIntegral . _unOb
+{-# INLINE unOb #-}
+
+
+-- | An atomic label.
+newtype Lb = Lb { _unLb :: Int16 }
+    deriving ( Show, Eq, Ord, Binary, A.IArray A.UArray
+             , G.Vector U.Vector, G.MVector U.MVector, U.Unbox
+             , Num, Ix, Bounds)
+
+
+-- | Smart label constructor.
+mkLb :: Int -> Lb
+mkLb = Lb . fromIntegral
+{-# INLINE mkLb #-}
+
+
+-- | Deconstract label.
+unLb :: Lb -> Int
+unLb = fromIntegral . _unLb
+{-# INLINE unLb #-}
+
+
+-- | An index of the label.
+type CbIx = Int
+
+
+-- | A feature index.  To every model feature a unique index is assigned.
+newtype FeatIx = FeatIx { _unFeatIx :: Int32 }
+    deriving ( Show, Eq, Ord, Binary, A.IArray A.UArray
+             , G.Vector U.Vector, G.MVector U.MVector, U.Unbox )
+
+
+-- | Smart feature index constructor.
+mkFeatIx :: Int -> FeatIx
+mkFeatIx = FeatIx . fromIntegral
+{-# INLINE mkFeatIx #-}
+
+
+-- | Deconstract feature index.
+unFeatIx :: FeatIx -> Int
+unFeatIx = fromIntegral . _unFeatIx
+{-# INLINE unFeatIx #-}
+
+
+----------------------------------------------------------------
+-- Complex label
+----------------------------------------------------------------
+
+
+-- TODO: Do we gain anything by representing the
+-- complex label with a byte array?  Complex labels
+-- should not be directly stored in a model, so if
+-- there is something to gain here, its not obvious.
+--
+-- Perhaps a list representation would be sufficient?
+
+
+-- -- | A complex label is an array of atomic labels.
+-- newtype Cb = Cb { unCb :: BA.ByteArray }
+
+
+-- | A complex label is a vector of atomic labels.
+newtype Cb = Cb { _unCb :: U.Vector Lb }
+    deriving (Show, Eq, Ord, Binary)
+
+
+-- | Smart complex label constructor.
+mkCb :: [Lb] -> Cb
+mkCb = Cb . U.fromList
+
+
+-- | Deconstract complex label.
+unCb :: Cb -> [Lb]
+unCb = U.toList . _unCb
+
+
+----------------------------------------------------------------
+-- Internal dataset representation
+----------------------------------------------------------------
+
+
+-- | A word is represented by a list of its observations
+-- and a list of its potential label interpretations.
+data X = X {
+    -- | A set of observations.
+      _unX :: U.Vector Ob
+    -- | A vector of potential labels.
+    , _unR :: V.Vector Cb }
+    deriving (Show, Eq, Ord)
+
+
+instance Binary X where
+    put X{..} = put _unX >> put _unR
+    get = X <$> get <*> get
+
+
+-- | Sentence of words.
+type Xs = V.Vector X
+
+
+-- | Smart `X` constructor.
+mkX :: [Ob] -> [Cb] -> X
+mkX x r = X (U.fromList x) (V.fromList r)
+{-# INLINE mkX #-}
+
+
+-- | List of observations.
+unX :: X -> [Ob]
+unX = U.toList . _unX
+{-# INLINE unX #-}
+
+
+-- | List of potential labels.
+unR :: X -> [Cb]
+unR = V.toList . _unR
+{-# INLINE unR #-}
+
+
+-- | Vector of chosen labels together with
+-- corresponding probabilities in log domain.
+newtype Y = Y { _unY :: V.Vector (Cb, Double) }
+    deriving (Show, Eq, Ord, Binary)
+
+
+-- | Y constructor.
+mkY :: [(Cb, Double)] -> Y
+mkY = Y . V.fromList . map (second log)
+{-# INLINE mkY #-}
+
+
+-- | Y deconstructor symetric to mkY.
+unY :: Y -> [(Cb, L.LogFloat)]
+unY = map (second L.logToLogFloat) . V.toList . _unY
+{-# INLINE unY #-}
+
+
+-- | Sentence of Y (label choices).
+type Ys = V.Vector Y
+
+
+-- | Potential label at the given position.
+lbAt :: X -> CbIx -> Cb
+lbAt x = (_unR x V.!)
+{-# INLINE lbAt #-}
+
+
+lbVec :: Xs -> Int -> V.Vector Cb
+lbVec xs = _unR . (xs V.!)
+{-# INLINE lbVec #-}
+
+
+-- | Number of potential labels at the given position of the sentence.
+lbNumI :: Xs -> Int -> Int
+lbNumI xs = V.length . lbVec xs
+{-# INLINE lbNumI #-}
+
+
+-- | Potential label at the given position and at the given index.
+lbOnI :: Xs -> Int -> CbIx -> Cb
+lbOnI xs = (V.!) . lbVec xs
+{-# INLINE lbOnI #-}
+
+
+-- | List of label indices at the given position.
+lbIxsI :: Xs -> Int -> [CbIx]
+lbIxsI xs i = [0 .. lbNum xs i - 1]
+{-# INLINE lbIxsI #-}
+
+
+-- | Number of potential labels at the given position of the sentence.
+-- Function extended to indices outside the positions' domain.
+lbNum :: Xs -> Int -> Int
+lbNum xs i
+    | i < 0 || i >= n   = 1
+    | otherwise         = lbNumI xs i
+  where
+    n = V.length xs
+{-# INLINE lbNum #-}
+
+
+-- | Potential label at the given position and at the given index.
+-- Return Nothing for positions outside the domain.
+lbOn :: Xs -> Int -> CbIx -> Maybe Cb
+lbOn xs i
+    | i < 0 || i >= n   = const Nothing
+    | otherwise         = Just . lbOnI xs i
+  where
+    n = V.length xs
+{-# INLINE lbOn #-}
+
+
+-- | List of label indices at the given position.  Function extended to
+-- indices outside the positions' domain.
+lbIxs :: Xs -> Int -> [CbIx]
+lbIxs xs i
+    | i < 0 || i >= n   = [0]
+    | otherwise         = lbIxsI xs i
+  where
+    n = V.length xs
+{-# INLINE lbIxs #-}
diff --git a/src/Data/CRF/Chain2/Tiers/Feature.hs b/src/Data/CRF/Chain2/Tiers/Feature.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Feature.hs
@@ -0,0 +1,206 @@
+module Data.CRF.Chain2.Tiers.Feature
+(
+-- * Feature
+  Feat (..)
+
+-- * Feature generation
+, obFeats
+, trFeats1
+, trFeats2
+, trFeats3
+
+-- * Feature extraction
+, presentFeats
+, hiddenFeats
+, obFeatsOn
+, trFeatsOn
+
+-- * Feature selection
+, FeatSel
+, selectPresent
+, selectHidden
+) where
+
+
+import           Control.Applicative ((<*>), (<$>))
+import           Data.Maybe (maybeToList)
+import           Data.List (zip4)
+import           Data.Binary (Binary, put, get, putWord8, getWord8)
+import qualified Data.Vector as V
+import qualified Data.Number.LogFloat as L
+
+import Data.CRF.Chain2.Tiers.Dataset.Internal
+
+
+----------------------------------------------------
+-- Feature
+----------------------------------------------------
+
+
+-- | Feature; every feature is associated to a layer with `ln` identifier.
+data Feat
+    -- | Second-order transition feature.
+    = TFeat3
+        { x1    :: {-# UNPACK #-} !Lb
+        , x2    :: {-# UNPACK #-} !Lb
+        , x3    :: {-# UNPACK #-} !Lb
+        , ln    :: {-# UNPACK #-} !Int }
+    -- | First-order transition feature.
+    | TFeat2
+        { x1    :: {-# UNPACK #-} !Lb
+        , x2    :: {-# UNPACK #-} !Lb
+        , ln    :: {-# UNPACK #-} !Int }
+    -- | Zero-order transition feature.
+    | TFeat1
+        { x1    :: {-# UNPACK #-} !Lb
+        , ln    :: {-# UNPACK #-} !Int }
+    -- | Observation feature.
+    | OFeat
+        { ob    :: {-# UNPACK #-} !Ob
+        , x1    :: {-# UNPACK #-} !Lb
+        , ln    :: {-# UNPACK #-} !Int }
+    deriving (Show, Eq, Ord)
+
+
+instance Binary Feat where
+    put (OFeat o x k)       = putWord8 0 >> put o >> put x >> put k
+    put (TFeat3 x y z k)    = putWord8 1 >> put x >> put y >> put z >> put k
+    put (TFeat2 x y k)      = putWord8 2 >> put x >> put y >> put k
+    put (TFeat1 x k)        = putWord8 3 >> put x >> put k
+    get = getWord8 >>= \i -> case i of
+        0   -> OFeat  <$> get <*> get <*> get
+        1   -> TFeat3 <$> get <*> get <*> get <*> get
+        2   -> TFeat2 <$> get <*> get <*> get
+        3   -> TFeat1 <$> get <*> get
+        _   -> error "get feature: unknown code"
+
+
+----------------------------------------------------
+-- Features generation
+----------------------------------------------------
+
+
+-- | Generate observation features.
+obFeats :: Ob -> Cb -> [Feat]
+obFeats ob' xs =
+    [ OFeat ob' x k
+    | (x, k) <- zip (unCb xs) [0..] ]
+
+
+-- | Generate zero-order transition features.
+trFeats1 :: Cb -> [Feat]
+trFeats1 xs =
+    [ TFeat1 x k
+    | (x, k) <- zip (unCb xs) [0..] ]
+
+
+-- | Generate first-order transition features.
+trFeats2 :: Cb -> Cb -> [Feat]
+trFeats2 xs1 xs2 =
+    [ TFeat2 x1' x2' k
+    | (x1', x2', k) <- zip3 (unCb xs1) (unCb xs2) [0..] ]
+
+
+-- | Generate second-order transition features.
+trFeats3 :: Cb -> Cb -> Cb -> [Feat]
+trFeats3 xs1 xs2 xs3 =
+    [ TFeat3 x1' x2' x3' k
+    | (x1', x2', x3', k) <- zip4 (unCb xs1) (unCb xs2) (unCb xs3) [0..] ]
+
+
+----------------------------------------------------
+-- Feature extraction
+----------------------------------------------------
+
+
+-- | Features present in the dataset element together with corresponding
+-- occurence probabilities.
+presentFeats :: Xs -> Ys -> [(Feat, L.LogFloat)]
+presentFeats xs ys = concat
+    [ obFs i ++ trFs i
+    | i <- [0 .. V.length xs - 1] ]
+  where
+    obFs i =
+        [ (ft, pr)
+        | o <- unX (xs V.! i)
+        , (u, pr) <- unY (ys V.! i)
+        , ft <- obFeats o u ]
+    trFs 0 =
+        [ (ft, pr)
+        | (u, pr) <- unY (ys V.! 0)
+        , ft <- trFeats1 u ]
+    trFs 1 =
+        [ (ft, pr1 * pr2)
+        | (u, pr1) <- unY (ys V.! 1)
+        , (v, pr2) <- unY (ys V.! 0)
+        , ft <- trFeats2 u v ]
+    trFs i =
+        [ (ft, pr1 * pr2 * pr3)
+        | (u, pr1) <- unY (ys V.! i)
+        , (v, pr2) <- unY (ys V.! (i-1))
+        , (w, pr3) <- unY (ys V.! (i-2))
+        , ft <- trFeats3 u v w ]
+
+
+-- | Features hidden in the dataset element.
+hiddenFeats :: Xs -> [Feat]
+hiddenFeats xs =
+    obFs ++ trFs
+  where
+    obFs = concat
+        [ obFeatsOn xs i u
+        | i <- [0 .. V.length xs - 1]
+        , u <- lbIxs xs i ]
+    trFs = concat
+        [ trFeatsOn xs i u v w
+        | i <- [0 .. V.length xs - 1]
+        , u <- lbIxs xs i
+        , v <- lbIxs xs $ i - 1
+        , w <- lbIxs xs $ i - 2 ]
+
+
+-- | Observation features on a given position and with respect
+-- to a given label (determined by idenex).
+obFeatsOn :: Xs -> Int -> CbIx -> [Feat]
+obFeatsOn xs i u = concat
+    [ obFeats ob' e
+    | e   <- lbs
+    , ob' <- unX (xs V.! i) ]
+  where 
+    lbs     = maybeToList (lbOn xs i u)
+{-# INLINE obFeatsOn #-}
+
+
+-- | Transition features on a given position and with respect
+-- to a given labels (determined by indexes).
+trFeatsOn :: Xs -> Int -> CbIx -> CbIx -> CbIx -> [Feat]
+trFeatsOn xs i u' v' w' =
+    doIt a b c
+  where
+    a = lbOn xs i       u'
+    b = lbOn xs (i - 1) v'
+    c = lbOn xs (i - 2) w'
+    doIt (Just u) (Just v) (Just w) = trFeats3 u v w
+    doIt (Just u) (Just v) _        = trFeats2 u v
+    doIt (Just u) _ _               = trFeats1 u
+    doIt _ _ _                      = []
+{-# INLINE trFeatsOn #-}
+
+
+----------------------------------------------------
+-- Feature selection
+----------------------------------------------------
+
+
+-- | A feature selection function type.
+type FeatSel = Xs -> Ys -> [Feat]
+
+
+-- | The 'presentFeats' adapted to fit feature selection specs.
+selectPresent :: FeatSel
+selectPresent xs = map fst . presentFeats xs
+
+
+-- | The 'hiddenFeats' adapted to fit feature selection specs.
+selectHidden :: FeatSel
+selectHidden xs _ = hiddenFeats xs
diff --git a/src/Data/CRF/Chain2/Tiers/Inference.hs b/src/Data/CRF/Chain2/Tiers/Inference.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Inference.hs
@@ -0,0 +1,229 @@
+{-# LANGUAGE RecordWildCards #-}
+
+module Data.CRF.Chain2.Tiers.Inference
+( tag
+, probs
+, marginals
+, expectedFeatures
+, accuracy
+, zx
+, zx'
+) where
+
+import           Data.Ord (comparing)
+import           Data.List (maximumBy)
+import qualified Data.Array as A
+import qualified Data.Vector as V
+import qualified Data.Number.LogFloat as L
+
+import           Control.Parallel.Strategies (rseq, parMap)
+import           Control.Parallel (par, pseq)
+import           GHC.Conc (numCapabilities)
+
+import           Data.CRF.Chain2.Tiers.Dataset.Internal
+import           Data.CRF.Chain2.Tiers.Model
+import           Data.CRF.Chain2.Tiers.Feature
+import           Data.CRF.Chain2.Tiers.Util (partition)
+import qualified Data.CRF.Chain2.Tiers.DP as DP
+
+
+-- Interface on top of internal implementation
+
+-- | Accumulation function.
+type AccF = [L.LogFloat] -> L.LogFloat
+
+type ProbArray = CbIx -> CbIx -> CbIx -> L.LogFloat
+
+computePsi :: Model -> Xs -> Int -> CbIx -> L.LogFloat
+computePsi crf xs i = (A.!) $ A.array (0, lbNum xs i - 1)
+    [ (k, onWord crf xs i k)
+    | k <- lbIxs xs i ]
+
+forward :: AccF -> Model -> Xs -> ProbArray
+forward acc crf sent = alpha where
+    alpha = DP.flexible3 (-1, V.length sent - 1)
+                (\i   -> (0, lbNum sent i - 1))
+                (\i _ -> (0, lbNum sent (i - 1) - 1))
+                (\t i -> withMem (computePsi crf sent i) t i)
+    withMem psi alpha i j k
+        | i == -1 = 1.0
+        | otherwise = acc
+            [ alpha (i - 1) k h * psi j
+            * onTransition crf sent i j k h
+            | h <- lbIxs sent (i - 2) ]
+
+backward :: AccF -> Model -> Xs -> ProbArray
+backward acc crf sent = beta where
+    beta = DP.flexible3 (0, V.length sent)
+               (\i   -> (0, lbNum sent (i - 1) - 1))
+               (\i _ -> (0, lbNum sent (i - 2) - 1))
+               (\t i -> withMem (computePsi crf sent i) t i)
+    withMem psi beta i j k
+        | i == V.length sent = 1.0
+        | otherwise = acc
+            [ beta (i + 1) h j * psi h
+            * onTransition crf sent i h j k
+            | h <- lbIxs sent i ]
+
+zxBeta :: ProbArray -> L.LogFloat
+zxBeta beta = beta 0 0 0
+
+zxAlpha :: AccF -> Xs -> ProbArray -> L.LogFloat
+zxAlpha acc sent alpha = acc
+    [ alpha (n - 1) i j
+    | i <- lbIxs sent (n - 1)
+    , j <- lbIxs sent (n - 2) ]
+    where n = V.length sent
+
+-- | Normalization factor computed for the Xs sentence using the
+-- backward computation.
+zx :: Model -> Xs -> L.LogFloat
+zx crf = zxBeta . backward sum crf
+
+-- | Normalization factor computed for the Xs sentence using the
+-- forward computation.
+zx' :: Model -> Xs -> L.LogFloat
+zx' crf sent = zxAlpha sum sent (forward sum crf sent)
+
+argmax :: (Ord b) => (a -> b) -> [a] -> (a, b)
+argmax f l = foldl1 choice $ map (\x -> (x, f x)) l
+    where choice (x1, v1) (x2, v2)
+              | v1 > v2 = (x1, v1)
+              | otherwise = (x2, v2)
+
+tagIxs :: Model -> Xs -> [Int]
+tagIxs crf sent = collectMaxArg (0, 0, 0) [] mem where
+    mem = DP.flexible3 (0, V.length sent)
+                       (\i   -> (0, lbNum sent (i - 1) - 1))
+                       (\i _ -> (0, lbNum sent (i - 2) - 1))
+                       (\t i -> withMem (computePsi crf sent i) t i)
+    withMem psiMem mem i j k
+        | i == V.length sent = (-1, 1)
+        | otherwise = argmax eval $ lbIxs sent i
+        where eval h =
+                  (snd $ mem (i + 1) h j) * psiMem h
+                  * onTransition crf sent i h j k
+    collectMaxArg (i, j, k) acc mem =
+        collect $ mem i j k
+        where collect (h, _)
+                  | h == -1 = reverse acc
+                  | otherwise = collectMaxArg (i + 1, h, j) (h:acc) mem
+
+-- | Find the most probable label sequence satisfying the constraints
+-- imposed over label values.
+tag :: Model -> Xs -> [Cb]
+tag crf sent =
+    let ixs = tagIxs crf sent
+    in  [lbAt x i | (x, i) <- zip (V.toList sent) ixs]
+
+-- | Tag potential labels with corresponding probabilities.
+probs :: Model -> Xs -> [[L.LogFloat]]
+probs crf sent =
+    let alpha = forward maximum crf sent
+        beta = backward maximum crf sent
+        normalize xs =
+            let d = - sum xs
+            in map (*d) xs
+        m1 k x = maximum
+            [ alpha k x y * beta (k + 1) x y
+            | y <- lbIxs sent (k - 1) ]
+    in  [ normalize [m1 i k | k <- lbIxs sent i]
+        | i <- [0 .. V.length sent - 1] ]
+
+-- | Tag potential labels with marginal probabilities.
+marginals :: Model -> Xs -> [[L.LogFloat]]
+marginals crf sent =
+    let alpha = forward sum crf sent
+        beta = backward sum crf sent
+    in  [ [ prob1 crf alpha beta sent i k
+          | k <- lbIxs sent i ]
+        | i <- [0 .. V.length sent - 1] ]
+
+goodAndBad :: Model -> Xs -> Ys -> (Int, Int)
+goodAndBad crf xs ys =
+    foldl gather (0, 0) $ zip labels labels'
+  where
+    labels  = [ (best . unY) (ys V.! i)
+              | i <- [0 .. V.length ys - 1] ]
+    best zs
+        | null zs   = Nothing
+        | otherwise = Just . fst $ maximumBy (comparing snd) zs
+    labels' = map Just $ tag crf xs
+    gather (good, bad) (x, y)
+        | x == y = (good + 1, bad)
+        | otherwise = (good, bad + 1)
+
+goodAndBad' :: Model -> [(Xs, Ys)] -> (Int, Int)
+goodAndBad' crf dataset =
+    let add (g, b) (g', b') = (g + g', b + b')
+    in  foldl add (0, 0) [goodAndBad crf x y | (x, y) <- dataset]
+
+-- | Compute the accuracy of the model with respect to the labeled dataset.
+accuracy :: Model -> [(Xs, Ys)] -> Double
+accuracy crf dataset =
+    let k = numCapabilities
+    	parts = partition k dataset
+        xs = parMap rseq (goodAndBad' crf) parts
+        (good, bad) = foldl add (0, 0) xs
+        add (g, b) (g', b') = (g + g', b + b')
+    in  fromIntegral good / fromIntegral (good + bad)
+
+prob3
+    :: Model -> ProbArray -> ProbArray -> Xs
+    -> Int -> (CbIx -> L.LogFloat) -> CbIx -> CbIx -> CbIx
+    -> L.LogFloat
+prob3 crf alpha beta sent k psiMem x y z =
+    alpha (k - 1) y z * beta (k + 1) x y * psiMem x
+    * onTransition crf sent k x y z / zxBeta beta
+{-# INLINE prob3 #-}
+
+prob2
+    :: Model -> ProbArray -> ProbArray
+    -> Xs -> Int -> CbIx -> CbIx -> L.LogFloat
+prob2 _ alpha beta _ k x y =
+    alpha k x y * beta (k + 1) x y / zxBeta beta
+{-# INLINE prob2 #-}
+
+prob1
+    :: Model -> ProbArray -> ProbArray
+    -> Xs -> Int -> CbIx -> L.LogFloat
+prob1 crf alpha beta sent k x = sum
+    [ prob2 crf alpha beta sent k x y
+    | y <- lbIxs sent (k - 1) ]
+
+expectedFeaturesOn
+    :: Model -> ProbArray -> ProbArray
+    -> Xs -> Int -> [(Feat, L.LogFloat)]
+expectedFeaturesOn crf alpha beta sent k =
+    fs3 ++ fs1
+    where psi = computePsi crf sent k
+          pr1 = prob1 crf alpha beta sent k
+          pr3 = prob3 crf alpha beta sent k psi
+          fs1 = [ (ft, pr) 
+                | a <- lbIxs sent k
+                , let pr = pr1 a
+                , ft <- obFs a ]
+    	  fs3 = [ (ft, pr) 
+                | a <- lbIxs sent k
+                , b <- lbIxs sent $ k - 1
+                , c <- lbIxs sent $ k - 2
+                , let pr = pr3 a b c
+                , ft <- trFs a b c ]
+          obFs = obFeatsOn sent k
+          trFs = trFeatsOn sent k
+
+-- | A list of features (represented by feature indices) defined within
+-- the context of the sentence accompanied by expected probabilities
+-- determined on the basis of the model.
+--
+-- One feature can occur multiple times in the output list.
+expectedFeatures :: Model -> Xs -> [(Feat, L.LogFloat)]
+expectedFeatures crf sent =
+    -- force parallel computation of alpha and beta tables
+    zx1 `par` zx2 `pseq` zx1 `pseq` concat
+      [ expectedFeaturesOn crf alpha beta sent k
+      | k <- [0 .. V.length sent - 1] ]
+    where alpha = forward sum crf sent
+          beta = backward sum crf sent
+          zx1 = zxAlpha sum sent alpha
+          zx2 = zxBeta beta
diff --git a/src/Data/CRF/Chain2/Tiers/Model.hs b/src/Data/CRF/Chain2/Tiers/Model.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Model.hs
@@ -0,0 +1,330 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE BangPatterns #-}
+
+
+module Data.CRF.Chain2.Tiers.Model
+( 
+-- * Model
+  Model (..)
+, mkModel
+, fromSet
+, fromMap
+, toMap
+
+-- * Potential
+, phi
+, index
+, onWord
+, onTransition
+) where
+
+
+import           Control.Applicative ((<$>), (<*>))
+import           Control.Monad (guard)
+import           Data.Binary (Binary, get, put)
+import           Data.Int (Int32)
+import qualified Data.Set as S
+import qualified Data.Map as M
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import           Data.Vector.Binary ()
+import qualified Data.Number.LogFloat as L
+
+import           Data.CRF.Chain2.Tiers.Dataset.Internal
+import           Data.CRF.Chain2.Tiers.Feature
+import qualified Data.CRF.Chain2.Tiers.Array as A
+
+
+-- | Dummy feature index.
+dummy :: FeatIx
+dummy = FeatIx (-1)
+{-# INLINE dummy #-}
+
+
+----------------------------------------------------
+-- (Feature -> FeatIx) map
+----------------------------------------------------
+
+
+-- | Transition map restricted to a particular tagging layer.
+-- Transitions of x form.
+type T1Map = A.Array Lb FeatIx
+
+
+-- | Transition map restricted to a particular tagging layer.
+-- Transitions of (x, y) form.
+type T2Map = A.Array (Lb, Lb) FeatIx
+
+
+-- | Transition map restricted to a particular tagging layer.
+-- Transitions of (x, y, z) form.
+type T3Map = A.Array (Lb, Lb, Lb) FeatIx
+
+
+mkT3Map :: [(Feat, FeatIx)] -> T3Map
+mkT3Map xs =
+    let ys = [((x, y, z), ix) | (TFeat3 x y z _, ix) <- xs]
+    in  A.mkArray dummy ys
+
+
+mkT2Map :: [(Feat, FeatIx)] -> T2Map
+mkT2Map xs =
+    let ys = [((x, y), ix) | (TFeat2 x y _, ix) <- xs]
+    in  A.mkArray dummy ys
+
+
+mkT1Map :: [(Feat, FeatIx)] -> T1Map
+mkT1Map xs =
+    let ys = [(x, ix) | (TFeat1 x _, ix) <- xs]
+    in  A.mkArray dummy ys
+
+
+unT3Map :: Int -> T3Map -> [(Feat, FeatIx)]
+unT3Map k t3 =
+    [ (TFeat3 x y z k, ix)
+    | ((x, y, z), ix) <- A.unArray t3
+    , ix /= dummy ]
+
+
+unT2Map :: Int -> T2Map -> [(Feat, FeatIx)]
+unT2Map k t2 =
+    [ (TFeat2 x y k, ix)
+    | ((x, y), ix) <- A.unArray t2
+    , ix /= dummy ]
+
+
+unT1Map :: Int -> T1Map -> [(Feat, FeatIx)]
+unT1Map k t1 =
+    [ (TFeat1 x k, ix)
+    | (x, ix) <- A.unArray t1
+    , ix /= dummy ]
+
+
+-- | Observation map restricted to a particular tagging layer.
+data OMap = OMap {
+    -- | Where memory ranges related to
+    -- individual observations begin?
+      oBeg  :: U.Vector Int32
+    -- | Labels related to individual observations.
+    , oLb   :: U.Vector Lb
+    -- | Feature indices related to individual (Ob, Lb) pairs.
+    , oIx   :: U.Vector FeatIx }
+
+
+instance Binary OMap where
+    put OMap{..} = put oBeg >> put oLb >> put oIx
+    get = OMap <$> get <*> get <*> get
+
+
+mkOMap :: [(Feat, FeatIx)] -> OMap
+mkOMap xs = OMap
+
+    { oBeg = U.fromList $ scanl (+) 0
+        [ fromIntegral (M.size lbMap)
+        | ob <- map mkOb [0 .. maxOb]
+        , let lbMap = maybe M.empty id $ M.lookup ob ftMap ]
+
+    , oLb = U.fromList . concat $
+        [ M.keys lbMap
+        | lbMap <- M.elems ftMap ]
+
+    , oIx = U.fromList . concat $
+        [ M.elems lbMap
+        | lbMap <- M.elems ftMap ] }
+
+  where
+
+    -- A feature map of type Map Ob (Map Lb FeatIx).
+    ftMap = fmap M.fromList $ M.fromListWith (++)
+        [ (ob, [(x, ix)])
+        | (OFeat ob x _, ix) <- xs ]
+
+    -- Max observation
+    maxOb = unOb . fst $ M.findMax ftMap
+
+
+-- | Deconstruct observation feature map given the layer identifier.
+unOMap :: Int -> OMap -> [(Feat, FeatIx)]
+unOMap k OMap{..} =
+    [ (OFeat o x k, i)
+    | (o, (p, q)) <- zip
+        (map mkOb [0..])
+        (pairs . map fromIntegral $ U.toList oBeg)
+    , (x, i)  <- zip
+        (U.toList $ U.slice p (q - p) oLb)
+        (U.toList $ U.slice p (q - p) oIx) ]
+  where
+    pairs xs = zip xs (tail xs)
+
+
+-- | Feature map restricted to a particular layer.
+data LayerMap = LayerMap
+    { t1Map     :: !T1Map
+    , t2Map     :: !T2Map
+    , t3Map     :: !T3Map
+    , obMap     :: !OMap }
+
+
+instance Binary LayerMap where
+    put LayerMap{..} = put t1Map >> put t2Map >> put t3Map >> put obMap
+    get = LayerMap <$> get <*> get <*> get <*> get
+
+
+-- | Deconstruct the layer map given the layer identifier.
+unLayerMap :: Int -> LayerMap -> [(Feat, FeatIx)]
+unLayerMap k LayerMap{..} = concat
+    [ unT1Map k t1Map
+    , unT2Map k t2Map
+    , unT3Map k t3Map
+    , unOMap  k obMap ]
+
+
+-- | Feature map is a vectro of layer maps.
+type FeatMap = V.Vector LayerMap
+
+
+-- | Get index of a feature.
+featIndex :: Feat -> FeatMap -> Maybe FeatIx
+featIndex (TFeat3 x y z k) v = do
+    m  <- t3Map <$> (v V.!? k)
+    ix <- m A.!? (x, y, z)
+    guard (ix /= dummy)
+    return ix
+featIndex (TFeat2 x y k) v = do
+    m  <- t2Map <$> (v V.!? k)
+    ix <- m A.!? (x, y)
+    guard (ix /= dummy)
+    return ix
+featIndex (TFeat1 x k) v = do
+    m  <- t1Map <$> (v V.!? k)
+    ix <- m A.!? x
+    guard (ix /= dummy)
+    return ix
+featIndex (OFeat ob x k) v = do
+    OMap{..} <- obMap <$> (v V.!? k)
+    p  <- fromIntegral <$> oBeg U.!? (unOb ob)
+    q  <- fromIntegral <$> oBeg U.!? (unOb ob + 1)
+    i  <- U.findIndex (==x) (U.slice p (q - p) oLb)
+    ix <- oIx U.!? (p + i)
+    -- guard (ix /= dummy)
+    return ix
+
+
+-- | Make feature map from a *set* of (feature, index) pairs.
+mkFeatMap :: [(Feat, FeatIx)] -> FeatMap
+mkFeatMap xs = V.fromList
+
+    -- TODO: We can first divide features between individual layers.
+    [ mkLayerMap $ filter (inLayer k . fst) xs
+    | k <- [0 .. maxLayerNum] ]
+
+  where
+
+    -- Make layer map.
+    mkLayerMap = LayerMap
+        <$> mkT1Map
+        <*> mkT2Map
+        <*> mkT3Map
+        <*> mkOMap
+
+    -- Number of layers (TODO: could be taken as mkFeatMap argument).
+    maxLayerNum = maximum $ map (ln.fst) xs
+
+    -- Check if feature is in a given layer.
+    inLayer k x | ln x == k     = True
+                | otherwise     = False
+
+
+-- | Deconstruct the feature map.
+unFeatMap :: FeatMap -> [(Feat, FeatIx)]
+unFeatMap fm = concat
+    [ unLayerMap i layer
+    | (i, layer) <- zip [0..] (V.toList fm) ]
+
+
+----------------------------------------------------
+-- Internal model
+----------------------------------------------------
+
+
+-- | Internal model data.
+data Model = Model
+    { values        :: U.Vector Double
+    , featMap       :: FeatMap }
+
+
+instance Binary Model where
+    put Model{..} = put values >> put featMap
+    get = Model <$> get <*> get
+
+
+-- | Construct model from a feature set.
+-- All values will be set to 1 in log domain.
+fromSet :: S.Set Feat -> Model
+fromSet ftSet = Model
+    { values    = U.replicate (S.size ftSet) 0.0
+    , featMap   =
+        let featIxs = map FeatIx [0..]
+            featLst = S.toList ftSet
+        in  mkFeatMap (zip featLst featIxs) }
+
+
+-- | Construct model from a dataset given a feature selection function.
+mkModel :: FeatSel -> [(Xs, Ys)] -> Model
+mkModel ftSel = fromSet . S.fromList . concatMap (uncurry ftSel)
+
+
+-- | Construct model from a (feature -> value) map.
+fromMap :: M.Map Feat L.LogFloat -> Model
+fromMap ftMap = Model
+    { values    = U.fromList . map L.logFromLogFloat $ M.elems ftMap
+    , featMap   =
+        let featIxs = map FeatIx [0..]
+            featLst = M.keys ftMap
+        in  mkFeatMap (zip featLst featIxs) }
+
+
+-- | Convert model to a (feature -> value) map.
+toMap :: Model -> M.Map Feat L.LogFloat
+toMap Model{..} = M.fromList
+    [ (ft, L.logToLogFloat (values U.! unFeatIx ix))
+    | (ft, ix) <- unFeatMap featMap ]
+
+
+----------------------------------------------------
+-- Potential
+----------------------------------------------------
+
+
+-- | Potential assigned to the feature -- exponential of the
+-- corresonding parameter.
+phi :: Model -> Feat -> L.LogFloat
+phi Model{..} ft = case featIndex ft featMap of
+    Just ix -> L.logToLogFloat (values U.! unFeatIx ix)
+    Nothing -> L.logToLogFloat (0 :: Float)
+{-# INLINE phi #-}
+
+
+-- | Index of a feature.
+index :: Model -> Feat -> Maybe FeatIx
+index Model{..} ft = featIndex ft featMap
+{-# INLINE index #-}
+
+
+-- | Observation potential on a given position and a
+-- given label (identified by index).
+onWord :: Model -> Xs -> Int -> CbIx -> L.LogFloat
+onWord crf xs i u =
+    product . map (phi crf) $ obFeatsOn xs i u
+{-# INLINE onWord #-}
+
+
+-- | Transition potential on a given position and a
+-- given labels (identified by indexes).
+onTransition :: Model -> Xs -> Int -> CbIx -> CbIx -> CbIx -> L.LogFloat
+onTransition crf xs i u w v =
+    product . map (phi crf) $ trFeatsOn xs i u w v
+{-# INLINE onTransition #-}
diff --git a/src/Data/CRF/Chain2/Tiers/Util.hs b/src/Data/CRF/Chain2/Tiers/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/CRF/Chain2/Tiers/Util.hs
@@ -0,0 +1,12 @@
+module Data.CRF.Chain2.Tiers.Util
+( partition
+) where
+
+import Data.List (transpose)
+
+partition :: Int -> [a] -> [[a]]
+partition n =
+    transpose . group n
+  where
+    group _ [] = []
+    group k xs = take k xs : (group k $ drop k xs)