diff --git a/Data/CRF/Chain1/Constrained.hs b/Data/CRF/Chain1/Constrained.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained.hs
@@ -0,0 +1,61 @@
+{-# LANGUAGE RecordWildCards #-}
+
+-- | The module provides first-order, linear-chain conditional random fields
+-- (CRFs) with position-wide constraints over label values.
+
+module Data.CRF.Chain1.Constrained
+(
+-- * Data types
+  Word (..)
+, unknown
+, Sent
+, Dist (unDist)
+, mkDist
+, WordL
+, annotate
+, SentL
+
+-- * CRF
+, CRF (..)
+-- ** Training
+, train
+-- ** Tagging
+, tag
+, tagK
+
+-- * Feature selection
+, hiddenFeats
+, presentFeats
+) where
+
+import Data.CRF.Chain1.Constrained.Dataset.External
+import Data.CRF.Chain1.Constrained.Dataset.Codec
+import Data.CRF.Chain1.Constrained.Feature.Present
+import Data.CRF.Chain1.Constrained.Feature.Hidden
+import Data.CRF.Chain1.Constrained.Train
+import qualified Data.CRF.Chain1.Constrained.Inference as I
+
+-- | Determine the most probable label sequence within the context of the
+-- given sentence using the model provided by the 'CRF'.
+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 ]
+
+-- | Determine the most probable label sets of the given size (at maximum)
+-- for each position in the input sentence.
+tagK :: (Ord a, Ord b) => Int -> CRF a b -> Sent a b -> [[b]]
+tagK k CRF{..} sent
+    = onWords . map decodeChoice
+    . I.tagK k model . encodeSent codec
+    $ sent
+  where
+    decodeChoice = decodeLabels codec . map fst
+    onWords xss =
+        [ take k $ unJusts codec word xs
+        | (word, xs) <- zip sent xss ]
diff --git a/Data/CRF/Chain1/Constrained/DP.hs b/Data/CRF/Chain1/Constrained/DP.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/DP.hs
@@ -0,0 +1,43 @@
+module Data.CRF.Chain1.Constrained.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/Data/CRF/Chain1/Constrained/Dataset/Codec.hs b/Data/CRF/Chain1/Constrained/Dataset/Codec.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Dataset/Codec.hs
@@ -0,0 +1,199 @@
+module Data.CRF.Chain1.Constrained.Dataset.Codec
+( Codec
+, CodecM
+
+, 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
+, unJusts
+) where
+
+import Control.Applicative ((<$>), (<*>), pure)
+import Data.Maybe (catMaybes, fromJust)
+import Data.Lens.Common (fstLens, sndLens)
+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.Chain1.Constrained.Dataset.Internal
+import Data.CRF.Chain1.Constrained.Dataset.External
+
+-- | A codec.  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, C.AtomCodec (Maybe b))
+
+-- | The empty codec.  The label part is initialized with Nothing
+-- member, which represents unknown labels.  It is taken on account
+-- in the model implementation because it is assigned to the
+-- lowest label code and the model assumes that the set of labels
+-- is of the {0, ..., 'lbMax'} form.
+empty :: Ord b => Codec a b
+empty =
+    let withNo = C.execCodec C.empty (C.encode C.idLens Nothing)
+    in  (C.empty, withNo)
+
+-- | Type synonym for the codec monad.  It is important to notice that by a
+-- codec we denote here a structure of two 'C.AtomCodec's while in the
+-- monad-codec package it denotes a 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 Ob . C.encode' fstLens
+
+-- | Encode the observation and do *not* update the codec.
+encodeObN :: Ord a => a -> CodecM a b (Maybe Ob)
+encodeObN = fmap (fmap Ob) . C.maybeEncode fstLens
+
+-- | Encode the label and update the codec.
+encodeLbU :: Ord b => b -> CodecM a b Lb
+encodeLbU = fmap Lb . C.encode sndLens . Just
+
+-- | Encode the label and do *not* update the codec.
+encodeLbN :: Ord b => b -> CodecM a b Lb
+encodeLbN x = do
+    my <- C.maybeEncode sndLens (Just x)
+    Lb <$> ( case my of
+        Just y  -> return y
+        Nothing -> fromJust <$> C.maybeEncode sndLens Nothing )
+
+-- | 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 . unDist) 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 . unDist) 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 = fromList . 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 the codec on the basis of the labeled dataset, return the
+-- resultant codec and the encoded dataset.
+mkCodec :: (Ord a, Ord b) => [SentL a b] -> (Codec a b, [(Xs, Ys)])
+mkCodec
+    = swap
+    . C.runCodec empty
+    . mapM encodeSentL'Cu
+  where
+    swap (x, y) = (y, x)
+
+-- | 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 codec = C.evalCodec codec . mapM encodeSentL'Cn
+
+-- | Encode the dataset with the codec.
+encodeData :: (Ord a, Ord b) => Codec a b -> [Sent a b] -> [Xs]
+encodeData codec = map (encodeSent codec)
+
+-- | Decode the label.
+decodeLabel :: Ord b => Codec a b -> Lb -> Maybe b
+decodeLabel codec x = C.evalCodec codec $ C.decode sndLens (unLb x)
+
+-- | Decode the sequence of labels.
+decodeLabels :: Ord b => Codec a b -> [Lb] -> [Maybe b]
+decodeLabels codec xs = C.evalCodec codec $
+    sequence [C.decode sndLens (unLb x) | x <- xs]
+
+hasLabel :: Ord b => Codec a b -> b -> Bool
+hasLabel codec x = M.member (Just x) (C.to $ snd codec)
+{-# INLINE hasLabel #-}
+
+-- | 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/Data/CRF/Chain1/Constrained/Dataset/External.hs b/Data/CRF/Chain1/Constrained/Dataset/External.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Dataset/External.hs
@@ -0,0 +1,62 @@
+module Data.CRF.Chain1.Constrained.Dataset.External
+( Word (..)
+, unknown
+, Sent
+, Dist (unDist)
+, mkDist
+, WordL
+, annotate
+, SentL
+) where
+
+import qualified Data.Set as S
+import qualified Data.Map as M
+
+-- | A Word is represented by a set of observations
+-- and a set of potential interpretation labels.
+-- When the set of potential labels is empty the word
+-- is considered to be unknown and the default potential
+-- set is used in its place.
+data Word a b = Word
+    { obs   :: S.Set a  -- ^ The set of observations
+    , lbs   :: S.Set b  -- ^ The set of potential interpretations.
+    } deriving (Show, Eq, Ord)
+
+-- | The word is considered to be unknown when the set of potential
+-- labels is empty.
+unknown :: Word a b -> Bool
+unknown word = S.size (lbs word) == 0
+{-# INLINE unknown #-}
+
+-- | 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 Dist a = Dist { unDist :: M.Map a Double }
+
+-- | Construct the probability distribution.
+mkDist :: Ord a => [(a, Double)] -> Dist a
+mkDist =
+    Dist . normalize . M.fromListWith (+)
+  where
+    normalize dist =
+        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, Dist b)
+
+-- | Annotate the word with the label.
+annotate :: Ord b => Word a b -> b -> WordL a b
+annotate w x
+    | x `S.member` lbs w    = (w, Dist (M.singleton x 1))
+    | otherwise             =
+        error "annotate: label not in the set of potential interpretations"
+
+-- | A sentence of labeled words.
+type SentL a b = [WordL a b]
diff --git a/Data/CRF/Chain1/Constrained/Dataset/Internal.hs b/Data/CRF/Chain1/Constrained/Dataset/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Dataset/Internal.hs
@@ -0,0 +1,109 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RecordWildCards #-}
+
+module Data.CRF.Chain1.Constrained.Dataset.Internal
+( Ob (..)
+, Lb (..)
+
+, X (..)
+, mkX
+, unX
+, unR
+, Xs
+
+, Y (..)
+, mkY
+, unY
+, Ys
+
+, AVec (unAVec)
+, fromList
+, fromSet
+) where
+
+import Data.Vector.Generic.Base
+import Data.Vector.Generic.Mutable
+import Data.Binary (Binary)
+import Data.Vector.Binary ()
+import Data.Ix (Ix)
+import qualified Data.Set as S
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+
+-- | An observation.
+newtype Ob = Ob { unOb :: Int }
+    deriving ( Show, Read, Eq, Ord, Binary
+             , Vector U.Vector, MVector U.MVector, U.Unbox )
+
+-- | A label.
+newtype Lb = Lb { unLb :: Int }
+    deriving ( Show, Read, Eq, Ord, Binary
+             , Vector U.Vector, MVector U.MVector, U.Unbox
+	     , Num, Ix )
+
+-- | Ascending vector of unique interger elements.
+newtype AVec a = AVec { unAVec :: U.Vector a }
+    deriving (Show, Read, Eq, Ord, Binary)
+
+-- | Smart AVec constructor which ensures that the
+-- underlying vector satisfies the AVec properties.
+fromList :: (Ord a, U.Unbox a) => [a] -> AVec a
+fromList = fromSet . S.fromList 
+{-# INLINE fromList #-}
+
+-- | Smart AVec constructor which ensures that the
+-- underlying vector satisfies the AVec properties.
+fromSet :: (Ord a, U.Unbox a) => S.Set a -> AVec a
+fromSet = AVec . U.fromList . S.toList 
+{-# INLINE fromSet #-}
+
+-- | A word represented by a list of its observations
+-- and a list of its potential label interpretations.
+data X
+    -- | The word with default set of potential interpretations.
+    = X { _unX :: AVec Ob }
+    -- | The word with custom set of potential labels.
+    | R { _unX :: AVec Ob
+        , _unR :: AVec Lb }
+    deriving (Show, Read, Eq, Ord)
+
+-- | X constructor.
+mkX :: [Ob] -> [Lb] -> X
+mkX x [] = X (fromList x)
+mkX x r  = R (fromList x) (fromList r)
+{-# INLINE mkX #-}
+
+-- | List of observations.
+unX :: X -> [Ob]
+unX = U.toList . unAVec . _unX
+{-# INLINE unX #-}
+
+-- | List of potential labels.
+unR :: AVec Lb -> X -> [Lb]
+unR r0 X{..} = U.toList . unAVec $ r0
+unR _  R{..} = U.toList . unAVec $ _unR
+{-# INLINE unR #-}
+
+-- | Sentence of words.
+type Xs = V.Vector X
+
+-- | Probability distribution over labels.  We assume, that when y is
+-- a member of chosen labels list it is also a member of the list
+-- potential labels for corresponding 'X' word.
+-- TODO: Perhaps we should substitute 'Lb's with label indices
+-- corresponding to labels from the vector of potential labels?
+newtype Y = Y { _unY :: AVec (Lb, Double) }
+    deriving (Show, Read, Eq, Ord)
+
+-- | Y constructor.
+mkY :: [(Lb, Double)] -> Y
+mkY = Y . fromList
+{-# INLINE mkY #-}
+
+-- | Y deconstructor symetric to mkY.
+unY :: Y -> [(Lb, Double)]
+unY = U.toList . unAVec . _unY
+{-# INLINE unY #-}
+
+-- | Sentence of Y (label choices).
+type Ys = V.Vector Y
diff --git a/Data/CRF/Chain1/Constrained/Feature.hs b/Data/CRF/Chain1/Constrained/Feature.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Feature.hs
@@ -0,0 +1,86 @@
+module Data.CRF.Chain1.Constrained.Feature
+( Feature (..)
+, isSFeat
+, isTFeat
+, isOFeat
+, featuresIn
+) where
+
+import Data.Binary (Binary, Get, put, get)
+import Control.Applicative ((<*>), (<$>))
+import qualified Data.Vector as V
+import qualified Data.Number.LogFloat as L
+
+import Data.CRF.Chain1.Constrained.Dataset.Internal
+
+-- | A Feature is either an observation feature OFeature o x, which
+-- models relation between observation o and label x assigned to
+-- the same word, or a transition feature TFeature x y (SFeature x
+-- for the first position in the sentence), which models relation
+-- between two subsequent labels, x (on i-th position) and y
+-- (on (i-1)-th positoin).
+data Feature
+    = SFeature
+        {-# UNPACK #-} !Lb
+    | TFeature
+        {-# UNPACK #-} !Lb
+        {-# UNPACK #-} !Lb
+    | OFeature
+        {-# UNPACK #-} !Ob
+        {-# UNPACK #-} !Lb
+    deriving (Show, Eq, Ord)
+
+instance Binary Feature where
+    put (SFeature x)   = put (0 :: Int) >> put x
+    put (TFeature x y) = put (1 :: Int) >> put (x, y)
+    put (OFeature o x) = put (2 :: Int) >> put (o, x)
+    get = do
+        k <- get :: Get Int
+        case k of
+            0 -> SFeature <$> get
+            1 -> TFeature <$> get <*> get
+            2 -> OFeature <$> get <*> get
+	    _ -> error "Binary Feature: unknown identifier"
+
+-- | Is it a 'SFeature'?
+isSFeat :: Feature -> Bool
+isSFeat (SFeature _) = True
+isSFeat _            = False
+{-# INLINE isSFeat #-}
+
+-- | Is it an 'OFeature'?
+isOFeat :: Feature -> Bool
+isOFeat (OFeature _ _) = True
+isOFeat _              = False
+{-# INLINE isOFeat #-}
+
+-- | Is it a 'TFeature'?
+isTFeat :: Feature -> Bool
+isTFeat (TFeature _ _) = True
+isTFeat _              = False
+{-# INLINE isTFeat #-}
+
+-- | Transition features with assigned probabilities for given position.
+trFeats :: Ys -> Int -> [(Feature, L.LogFloat)]
+trFeats ys 0 =
+    [ (SFeature x, L.logFloat px)
+    | (x, px) <- unY (ys V.! 0) ]
+trFeats ys k =
+    [ (TFeature x y, L.logFloat px * L.logFloat py)
+    | (x, px) <- unY (ys V.! k)
+    , (y, py) <- unY (ys V.! (k-1)) ]
+
+-- | Observation features with assigned probabilities for a given position.
+obFeats :: Xs -> Ys -> Int -> [(Feature, L.LogFloat)]
+obFeats xs ys k =
+    [ (OFeature o x, L.logFloat px)
+    | (x, px) <- unY (ys V.! k)
+    , o       <- unX (xs V.! k) ]
+
+-- | All features with assigned probabilities for given position.
+features :: Xs -> Ys -> Int -> [(Feature, L.LogFloat)]
+features xs ys k = trFeats ys k ++ obFeats xs ys k
+
+-- | All features with assigned probabilities in given labeled sentence.
+featuresIn :: Xs -> Ys -> [(Feature, L.LogFloat)]
+featuresIn xs ys = concatMap (features xs ys) [0 .. V.length xs - 1]
diff --git a/Data/CRF/Chain1/Constrained/Feature/Hidden.hs b/Data/CRF/Chain1/Constrained/Feature/Hidden.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Feature/Hidden.hs
@@ -0,0 +1,59 @@
+-- | The module provides feature selection functions which extract
+-- hidden features, i.e. all features which can be constructed 
+-- on the basis of observations and potential labels (constraints)
+-- corresponding to individual words.
+--
+-- You can mix functions defined here with the selection functions
+-- from the "Data.CRF.Chain1.Constrained.Feature.Present" module.
+
+module Data.CRF.Chain1.Constrained.Feature.Hidden
+( hiddenFeats
+, hiddenOFeats
+, hiddenTFeats
+, hiddenSFeats
+) where
+
+import qualified Data.Vector as V
+
+import Data.CRF.Chain1.Constrained.Dataset.Internal
+import Data.CRF.Chain1.Constrained.Feature
+
+-- | Hidden 'OFeature's which can be constructed based on the dataset.
+-- The default set of potential interpretations is used for all unknown words.
+hiddenOFeats :: AVec Lb -> [(Xs, b)] -> [Feature]
+hiddenOFeats r0 ds =
+    concatMap f ds
+  where
+    f = concatMap oFeats . V.toList . fst
+    oFeats x =
+        [ OFeature o y
+        | o <- unX x
+        , y <- unR r0 x ]
+
+-- | Hidden 'TFeature's which can be constructed based on the dataset.
+-- The default set of potential interpretations is used for all unknown words.
+hiddenTFeats :: AVec Lb -> [(Xs, b)] -> [Feature]
+hiddenTFeats r0 ds =
+    concatMap (tFeats . fst) ds
+  where
+    tFeats xs = concatMap (tFeatsOn xs) [1 .. V.length xs - 1]
+    tFeatsOn xs k =
+        [ TFeature x y
+        | x <- unR r0 (xs V.! k)
+        , y <- unR r0 (xs V.! (k-1)) ]
+
+-- | Hidden 'SFeature's which can be constructed based on the dataset.
+-- The default set of potential interpretations is used for all unknown words.
+hiddenSFeats :: AVec Lb -> [(Xs, b)] -> [Feature]
+hiddenSFeats r0 ds =
+    let sFeats xs = [SFeature x | x <- unR r0 (xs V.! 0)]
+    in  concatMap (sFeats . fst) ds
+
+-- | Hidden 'Feature's of all types which can be constructed
+-- on the basis of the dataset.  The default set of potential
+-- interpretations is used for all unknown words.
+hiddenFeats :: AVec Lb -> [(Xs, b)] -> [Feature]
+hiddenFeats r0 ds
+    =  hiddenOFeats r0 ds
+    ++ hiddenTFeats r0 ds
+    ++ hiddenSFeats r0 ds
diff --git a/Data/CRF/Chain1/Constrained/Feature/Present.hs b/Data/CRF/Chain1/Constrained/Feature/Present.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Feature/Present.hs
@@ -0,0 +1,56 @@
+-- | The module provides feature selection functions which extract
+-- features present in the dataset, i.e. features which directly occure
+-- the dataset.
+--
+-- You can mix functions defined here with the selection functions
+-- from the "Data.CRF.Chain1.Constrained.Feature.Hidden" module.
+
+module Data.CRF.Chain1.Constrained.Feature.Present
+( presentFeats
+, presentOFeats
+, presentTFeats
+, presentSFeats
+) where
+
+import qualified Data.Vector as V
+
+import Data.CRF.Chain1.Constrained.Dataset.Internal
+import Data.CRF.Chain1.Constrained.Feature
+
+-- | 'OFeature's which occur in the dataset.
+presentOFeats :: [(Xs, Ys)] -> [Feature]
+presentOFeats ds =
+    concatMap sentOFeats ds
+  where
+    sentOFeats (xs, ys) = concatMap oFeatsOn (zip (V.toList xs) (V.toList ys))
+    oFeatsOn (x, choice) =
+        [ OFeature o y
+        | o <- unX x
+        , y <- lbs choice ]
+
+-- | 'TFeature's which occur in the dataset.
+presentTFeats :: [(a, Ys)] -> [Feature]
+presentTFeats ds =
+    concatMap (sentTFeats.snd) ds
+  where
+    sentTFeats ys = concatMap (tFeatsOn ys) [1 .. V.length ys - 1]
+    tFeatsOn ys k =
+        [ TFeature x y
+        | x <- lbs (ys V.! k)
+        , y <- lbs (ys V.! (k-1)) ]
+
+-- | 'SFeature's which occur in the dataset.
+presentSFeats :: [(a, Ys)] -> [Feature]
+presentSFeats ds =
+    let sentSFeats s = [SFeature x | x <- lbs (s V.! 0)] 
+    in  concatMap (sentSFeats.snd) ds
+
+-- | 'Feature's of all kinds which occur in the dataset.
+presentFeats :: [(Xs, Ys)] -> [Feature]
+presentFeats ds
+    =  presentOFeats ds
+    ++ presentTFeats ds
+    ++ presentSFeats ds
+
+lbs :: Y -> [Lb]
+lbs = map fst . unY
diff --git a/Data/CRF/Chain1/Constrained/Inference.hs b/Data/CRF/Chain1/Constrained/Inference.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Inference.hs
@@ -0,0 +1,247 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TupleSections #-}
+
+-- Inference with CRFs.
+
+module Data.CRF.Chain1.Constrained.Inference
+( tag
+, tagK
+, marginals
+, accuracy
+, expectedFeaturesIn
+, zx
+, zx'
+) where
+
+import Control.Applicative ((<$>))
+import Data.Maybe (catMaybes)
+import Data.List (maximumBy, sortBy)
+import Data.Function (on)
+import qualified Data.Array as A
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+
+import Control.Parallel.Strategies (rseq, parMap)
+import Control.Parallel (par, pseq)
+import GHC.Conc (numCapabilities)
+import qualified Data.Number.LogFloat as L
+
+import qualified Data.CRF.Chain1.Constrained.DP as DP
+import Data.CRF.Chain1.Constrained.Util (partition)
+import Data.CRF.Chain1.Constrained.Dataset.Internal
+import Data.CRF.Chain1.Constrained.Feature (Feature(..))
+import Data.CRF.Chain1.Constrained.Model
+import Data.CRF.Chain1.Constrained.Intersect
+
+type LbIx       = Int
+type ProbArray  = Int -> LbIx -> L.LogFloat
+
+-- Some basic definitions.
+
+-- | Vector of potential labels on the given position of the sentence.
+lbVec :: Model -> Xs -> Int -> AVec Lb
+lbVec crf xs k = case xs V.! k of
+    X _     -> (r0 crf)
+    R _ r   -> r
+{-# INLINE lbVec #-}
+
+-- | Number of potential labels on the given position of the sentence.
+lbNum :: Model -> Xs -> Int -> Int
+lbNum crf xs = (U.length . unAVec) . lbVec crf xs
+{-# INLINE lbNum #-}
+
+-- | Potential label on the given vector position.
+lbOn :: Model -> X -> Int -> Lb
+lbOn crf (X _)   = (unAVec (r0 crf) U.!)
+lbOn _   (R _ r) = (unAVec r U.!)
+{-# INLINE lbOn #-}
+
+lbIxs :: Model -> Xs -> Int -> [(Int, Lb)]
+lbIxs crf xs = zip [0..] . U.toList . unAVec . lbVec crf xs
+{-# INLINE lbIxs #-}
+
+-- | Compute the table of potential products associated with 
+-- observation features for the given sentence position.
+computePsi :: Model -> Xs -> Int -> LbIx -> L.LogFloat
+computePsi crf xs i = (A.!) $ A.accumArray (*) 1 bounds
+    [ (k, valueL crf ix)
+    | ob <- unX (xs V.! i)
+    , (k, ix) <- intersect (obIxs crf ob) (lbVec crf xs i) ]
+  where
+    bounds = (0, lbNum crf xs i - 1)
+
+-- | Forward table computation.
+forward :: Model -> Xs -> ProbArray
+forward crf xs = alpha where
+    alpha = DP.flexible2 (0, V.length xs) bounds
+        (\t i -> withMem (computePsi crf xs i) t i)
+    bounds i
+        | i == V.length xs = (0, 0)
+        | otherwise = (0, lbNum crf xs i - 1)
+    withMem psi alpha i
+        | i == V.length xs = const u
+        | i == 0 = \j ->
+            let x = lbOn crf (xs V.! i) j
+            in  psi j * sgValue crf x
+        | otherwise = \j ->
+            let x = lbOn crf (xs V.! i) j
+            in  psi j * ((u - v x) + w x)
+      where
+        u = sum
+            [ alpha (i-1) k
+            | (k, _) <- lbIxs crf xs (i-1) ]
+        v x = sum
+            [ alpha (i-1) k
+            | (k, _) <- intersect (prevIxs crf x) (lbVec crf xs (i-1)) ]
+        w x = sum
+            [ alpha (i-1) k * valueL crf ix
+            | (k, ix) <- intersect (prevIxs crf x) (lbVec crf xs (i-1)) ]
+
+-- | Backward table computation.
+backward :: Model -> Xs -> ProbArray
+backward crf xs = beta where
+    beta = DP.flexible2 (0, V.length xs) bounds
+        (\t i -> withMem (computePsi crf xs i) t i)
+    bounds i
+        | i == 0    = (0, 0)
+        | otherwise = (0, lbNum crf xs (i-1) - 1)
+    withMem psi beta i
+        | i == V.length xs = const 0
+        | i == 0 = const $ sum
+            [ beta (i+1) k * psi k
+            * sgValue crf (lbOn crf (xs V.! i) k)
+            | (k, _) <- lbIxs crf xs i ]
+        | otherwise = \j ->
+            let y = lbOn crf (xs V.! (i-1)) j
+            in  (u - v y) + w y
+      where
+        u = sum
+            [ beta (i+1) k * psi k
+            | (k, _ ) <- lbIxs crf xs i ]
+        v y = sum
+            [ beta (i+1) k * psi k
+            | (k, _ ) <- intersect (nextIxs crf y) (lbVec crf xs i) ]
+        w y = sum
+            [ beta (i+1) k * psi k * valueL crf ix
+            | (k, ix) <- intersect (nextIxs crf y) (lbVec crf xs i) ]
+
+zxBeta :: ProbArray -> L.LogFloat
+zxBeta beta = beta 0 0
+
+zxAlpha :: Xs -> ProbArray -> L.LogFloat
+zxAlpha xs alpha = alpha (V.length xs) 0
+
+-- | Normalization factor computed for the 'Xs' sentence using the
+-- backward computation.
+zx :: Model -> Xs -> L.LogFloat
+zx crf = zxBeta . backward crf
+
+-- | Normalization factor computed for the 'Xs' sentence using the
+-- forward computation.
+zx' :: Model -> Xs -> L.LogFloat
+zx' crf sent = zxAlpha sent (forward crf sent)
+
+-- | Tag probabilities with respect to marginal distributions.
+marginals :: Model -> Xs -> [[(Lb, L.LogFloat)]]
+marginals crf xs =
+    let alpha = forward crf xs
+        beta = backward crf xs
+    in  [ [ (x, prob1 alpha beta i k)
+          | (k, x) <- lbIxs crf xs i ]
+        | i <- [0 .. V.length xs - 1] ]
+
+-- | Get (at most) k best tags for each word and return them in
+-- descending order.  TODO: Tagging with respect to marginal
+-- distributions might not be the best idea.  Think of some
+-- more elegant method.
+tagK :: Int -> Model -> Xs -> [[(Lb, L.LogFloat)]]
+tagK k crf xs = map
+    ( take k
+    . reverse
+    . sortBy (compare `on` snd)
+    ) (marginals crf xs)
+
+-- | Find the most probable label sequence (with probabilities of individual
+-- lables determined with respect to marginal distributions) satisfying the
+-- constraints imposed over label values.
+tag :: Model -> Xs -> [Lb]
+tag crf = map (fst . head) . (tagK 1 crf)
+
+prob1 :: ProbArray -> ProbArray -> Int -> LbIx -> L.LogFloat
+prob1 alpha beta k x =
+    alpha k x * beta (k + 1) x / zxBeta beta
+{-# INLINE prob1 #-}
+
+prob2 :: Model -> ProbArray -> ProbArray -> Int -> (LbIx -> L.LogFloat)
+      -> LbIx -> LbIx -> FeatIx -> L.LogFloat
+prob2 crf alpha beta k psi x y ix
+    = alpha (k - 1) y * beta (k + 1) x
+    * psi x * valueL crf ix / zxBeta beta
+{-# INLINE prob2 #-}
+
+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 (compare `on` 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)
+
+expectedFeaturesOn
+    :: Model -> ProbArray -> ProbArray -> Xs
+    -> Int -> [(FeatIx, L.LogFloat)]
+expectedFeaturesOn crf alpha beta xs i =
+    tFeats ++ oFeats
+  where
+    psi = computePsi crf xs i
+    pr1 = prob1     alpha beta i
+    pr2 = prob2 crf alpha beta i psi
+
+    oFeats = [ (ix, pr1 k) 
+             | o <- unX (xs V.! i)
+             , (k, ix) <- intersect (obIxs crf o) (lbVec crf xs i) ]
+
+    tFeats
+        | i == 0 = catMaybes
+            [ (, pr1 k) <$> featToIx crf (SFeature x)
+            | (k, x) <- lbIxs crf xs i ]
+        | otherwise =
+            [ (ix, pr2 k l ix)
+            | (k,  x) <- lbIxs crf xs i
+            , (l, ix) <- intersect (prevIxs crf x) (lbVec crf xs (i-1)) ]
+
+-- | 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.
+expectedFeaturesIn :: Model -> Xs -> [(FeatIx, L.LogFloat)]
+expectedFeaturesIn crf xs = zxF `par` zxB `pseq` zxF `pseq`
+    concat [expectedOn k | k <- [0 .. V.length xs - 1] ]
+  where
+    expectedOn = expectedFeaturesOn crf alpha beta xs
+    alpha = forward crf xs
+    beta = backward crf xs
+    zxF = zxAlpha xs alpha
+    zxB = zxBeta beta
diff --git a/Data/CRF/Chain1/Constrained/Intersect.hs b/Data/CRF/Chain1/Constrained/Intersect.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Intersect.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE TupleSections #-}
+{-# LANGUAGE BangPatterns #-}
+
+module Data.CRF.Chain1.Constrained.Intersect
+( intersect
+) where
+
+import qualified Data.Vector.Unboxed as U
+
+import Data.CRF.Chain1.Constrained.Dataset.Internal (Lb, AVec, unAVec)
+import Data.CRF.Chain1.Constrained.Model (FeatIx)
+
+-- | Assumption: both input list are given in an ascending order.
+intersect
+    :: AVec (Lb, FeatIx)    -- ^ Vector of (label, features index) pairs
+    -> AVec Lb              -- ^ Vector of labels
+    -- | Intersection of arguments: vector indices from the second list
+    -- and feature indices from the first list.
+    -> [(Int, FeatIx)]
+intersect xs' ys'
+    | n == 0 || m == 0 = []
+    | otherwise = merge xs ys
+  where
+    xs = unAVec xs'
+    ys = unAVec ys'
+    n = U.length ys
+    m = U.length xs
+
+merge :: U.Vector (Lb, FeatIx) -> U.Vector Lb -> [(Int, FeatIx)]
+merge xs ys = doIt 0 0
+  where
+    m = U.length xs
+    n = U.length ys
+    doIt i j
+        | i >= m || j >= n = []
+        | otherwise = case compare x y of
+            EQ -> (j, ix) : doIt (i+1) (j+1)
+            LT -> doIt (i+1) j
+            GT -> doIt i (j+1)
+      where
+        (x, ix) = xs `U.unsafeIndex` i
+        y = ys `U.unsafeIndex` j
diff --git a/Data/CRF/Chain1/Constrained/Model.hs b/Data/CRF/Chain1/Constrained/Model.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Model.hs
@@ -0,0 +1,231 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+
+-- | Internal implementation of the CRF model.
+
+module Data.CRF.Chain1.Constrained.Model
+( FeatIx (..)
+, Model (..)
+, mkModel
+, valueL
+, featToIx
+, featToJustIx
+, featToJustInt
+, sgValue
+, sgIxs
+, obIxs
+, nextIxs
+, prevIxs
+) where
+
+import Control.Applicative ((<$>), (<*>))
+import Data.Maybe (fromJust)
+import Data.List (groupBy, sort)
+import Data.Function (on)
+import Data.Binary
+import qualified Data.Vector.Generic.Base as G
+import qualified Data.Vector.Generic.Mutable as G
+import qualified Data.Set as Set
+import qualified Data.Map as M
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector as V
+import qualified Data.Number.LogFloat as L
+
+import Data.CRF.Chain1.Constrained.Feature
+import Data.CRF.Chain1.Constrained.Dataset.Internal hiding (fromList)
+import qualified Data.CRF.Chain1.Constrained.Dataset.Internal as A
+
+-- | A feature index.  To every model feature a unique index is assigned.
+newtype FeatIx = FeatIx { unFeatIx :: Int }
+    deriving ( Show, Eq, Ord, Binary
+             , G.Vector U.Vector, G.MVector U.MVector, U.Unbox )
+
+-- | A label and a feature index determined by that label.
+type LbIx   = (Lb, FeatIx)
+
+dummyFeatIx :: FeatIx
+dummyFeatIx = FeatIx (-1)
+{-# INLINE dummyFeatIx #-}
+
+isDummy :: FeatIx -> Bool
+isDummy (FeatIx ix) = ix < 0
+{-# INLINE isDummy #-}
+
+notDummy :: FeatIx -> Bool
+notDummy = not . isDummy
+{-# INLINE notDummy #-}
+
+-- | The model is realy a map from features to potentials, but for the sake
+-- of efficiency the internal representation is more complex.
+data Model = Model {
+    -- | Value (potential) of the model for feature index.
+      values    :: U.Vector Double
+    -- | A map from features to feature indices
+    , ixMap     :: M.Map Feature FeatIx
+    -- | Default set of potential labels.
+    , r0        :: AVec Lb
+    -- | Singular feature index for the given label.  Index is equall to -1
+    -- if feature is not present in the model.
+    , sgIxsV 	:: U.Vector FeatIx
+    -- | Set of labels for the given observation which, together with the
+    -- observation, constitute an observation feature of the model. 
+    , obIxsV    :: V.Vector (AVec LbIx)
+    -- | Set of ,,previous'' labels for the value of the ,,current'' label.
+    -- Both labels constitute a transition feature present in the the model.
+    , prevIxsV  :: V.Vector (AVec LbIx)
+    -- | Set of ,,next'' labels for the value of the ,,current'' label.
+    -- Both labels constitute a transition feature present in the the model.
+    , nextIxsV  :: V.Vector (AVec LbIx) }
+
+instance Binary Model where
+    put crf = do
+        put $ values crf
+        put $ ixMap crf
+        put $ r0 crf
+        put $ sgIxsV crf
+        put $ obIxsV crf
+        put $ prevIxsV crf
+        put $ nextIxsV crf
+    get = Model <$> get <*> get <*> get <*> get <*> get <*> get <*> get
+
+-- | Construct CRF model from the associations list.  We assume that
+-- the set of labels is of the {0, 1, .. 'lbMax'} form and, similarly,
+-- the set of observations is of the {0, 1, .. 'obMax'} form.
+-- There should be no repetition of features in the input list.
+-- TODO: We can change this function to take M.Map Feature Double.
+fromList :: [(Feature, Double)] -> Model
+fromList fs =
+    let _ixMap = M.fromList $ zip
+            (map fst fs)
+            (map FeatIx [0..])
+    
+        sFeats = [feat | (feat, _val) <- fs, isSFeat feat]
+        tFeats = [feat | (feat, _val) <- fs, isTFeat feat]
+        oFeats = [feat | (feat, _val) <- fs, isOFeat feat]
+
+        obMax = (unOb . maximum . Set.toList . obSet) (map fst fs)
+        lbs   = (Set.toList . lbSet) (map fst fs)
+        lbMax = (unLb . maximum) lbs
+        _r0   = A.fromList lbs
+        
+        _sgIxsV = sgVects lbMax
+            [ (unLb x, featToJustIx crf feat)
+            | feat@(SFeature x) <- sFeats ]
+
+        _prevIxsV = adjVects lbMax
+            [ (unLb x, (y, featToJustIx crf feat))
+            | feat@(TFeature x y) <- tFeats ]
+
+        _nextIxsV = adjVects lbMax
+            [ (unLb y, (x, featToJustIx crf feat))
+            | feat@(TFeature x y) <- tFeats ]
+
+        _obIxsV = adjVects obMax
+            [ (unOb o, (x, featToJustIx crf feat))
+            | feat@(OFeature o x) <- oFeats ]
+
+        -- | Adjacency vectors.
+        adjVects n xs =
+            V.replicate n (A.fromList []) V.// update
+          where
+            update = map mkVect $ groupBy ((==) `on` fst) $ sort xs
+            mkVect (y:ys) = (fst y, A.fromList $ map snd (y:ys))
+            mkVect [] = error "mkVect: null list"
+
+        sgVects n xs = U.replicate n dummyFeatIx U.// xs
+
+        _values = U.replicate (length fs) 0.0
+            U.// [ (featToJustInt crf feat, val)
+                 | (feat, val) <- fs ]
+        crf = Model _values _ixMap _r0 _sgIxsV _obIxsV _prevIxsV _nextIxsV
+    in  crf
+
+-- | Compute the set of observations.
+obSet :: [Feature] -> Set.Set Ob
+obSet =
+    Set.fromList . concatMap toObs
+  where
+    toObs (OFeature o _) = [o]
+    toObs _              = []
+
+-- | Compute the set of labels.
+lbSet :: [Feature] -> Set.Set Lb
+lbSet =
+    Set.fromList . concatMap toLbs
+  where
+    toLbs (SFeature x)   = [x]
+    toLbs (OFeature _ x) = [x]
+    toLbs (TFeature x y) = [x, y]
+
+-- | Construct the model from the list of features.  All parameters will be
+-- set to 0.  There can be repetitions in the input list.
+-- We assume that the set of labels is of the {0, 1, .. 'lbMax'} form and,
+-- similarly, the set of observations is of the {0, 1, .. 'obMax'} form.
+mkModel :: [Feature] -> Model
+mkModel fs =
+    let fSet = Set.fromList fs
+        fs'  = Set.toList fSet
+        vs   = replicate (Set.size fSet) 0.0
+    in  fromList (zip fs' vs)
+
+-- | Model potential defined for the given feature interpreted as a
+-- number in logarithmic domain.
+valueL :: Model -> FeatIx -> L.LogFloat
+valueL crf (FeatIx i) = L.logToLogFloat (values crf U.! i)
+{-# INLINE valueL #-}
+
+-- | Determine index for the given feature.
+featToIx :: Model -> Feature -> Maybe FeatIx
+featToIx crf feat = M.lookup feat (ixMap crf)
+{-# INLINE featToIx #-}
+
+-- | Determine index for the given feature.  Throw error when
+-- the feature is not a member of the model. 
+featToJustIx :: Model -> Feature -> FeatIx
+featToJustIx _crf = fromJust . featToIx _crf
+{-# INLINE featToJustIx #-}
+
+-- | Determine index for the given feature and return it as an integer.
+-- Throw error when the feature is not a member of the model.
+featToJustInt :: Model -> Feature -> Int
+featToJustInt _crf = unFeatIx . featToJustIx _crf
+{-# INLINE featToJustInt #-}
+
+-- | Potential value (in log domain) of the singular feature with the
+-- given label.  The value defaults to 1 (0 in log domain) when the feature
+-- is not a member of the model.
+sgValue :: Model -> Lb -> L.LogFloat
+sgValue crf (Lb x) = 
+    case unFeatIx (sgIxsV crf U.! x) of
+        -- TODO: Is the value correct?
+        -1 -> L.logToLogFloat (0 :: Float)
+        ix -> L.logToLogFloat (values crf U.! ix)
+
+-- | List of labels which can be located on the first position of
+-- a sentence together with feature indices determined by them.
+sgIxs :: Model -> [LbIx]
+sgIxs crf
+    = filter (notDummy . snd)
+    . zip (map Lb [0..])
+    . U.toList $ sgIxsV crf
+{-# INLINE sgIxs #-}
+
+-- | List of labels which constitute a valid feature in combination with
+-- the given observation accompanied by feature indices determined by
+-- these labels.
+obIxs :: Model -> Ob -> AVec LbIx
+obIxs crf x = obIxsV crf V.! unOb x
+{-# INLINE obIxs #-}
+
+-- | List of ,,next'' labels which constitute a valid feature in combination
+-- with the ,,current'' label accompanied by feature indices determined by
+-- ,,next'' labels.
+nextIxs :: Model -> Lb -> AVec LbIx
+nextIxs crf x = nextIxsV crf V.! unLb x
+{-# INLINE nextIxs #-}
+
+-- | List of ,,previous'' labels which constitute a valid feature in
+-- combination with the ,,current'' label accompanied by feature indices
+-- determined by ,,previous'' labels.
+prevIxs :: Model -> Lb -> AVec LbIx
+prevIxs crf x = prevIxsV crf V.! unLb x
+{-# INLINE prevIxs #-}
diff --git a/Data/CRF/Chain1/Constrained/Train.hs b/Data/CRF/Chain1/Constrained/Train.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Train.hs
@@ -0,0 +1,108 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE PatternGuards #-}
+
+module Data.CRF.Chain1.Constrained.Train
+( CRF (..)
+, train
+) where
+
+import Control.Applicative ((<$>), (<*>))
+import System.IO (hSetBuffering, stdout, BufferMode (..))
+import Data.Binary (Binary, put, get)
+import qualified Data.Set as S
+import qualified Data.Map as M
+import qualified Data.Vector as V
+import qualified Numeric.SGD as SGD
+import qualified Numeric.SGD.LogSigned as L
+
+import Data.CRF.Chain1.Constrained.Dataset.Internal
+import Data.CRF.Chain1.Constrained.Dataset.External (SentL, unknown, unDist)
+import Data.CRF.Chain1.Constrained.Dataset.Codec
+    (mkCodec, Codec, encodeDataL, encodeLabels)
+import Data.CRF.Chain1.Constrained.Feature (Feature, featuresIn)
+import Data.CRF.Chain1.Constrained.Model
+    (Model (..), mkModel, FeatIx (..), featToJustInt)
+import Data.CRF.Chain1.Constrained.Inference (accuracy, expectedFeaturesIn)
+
+-- | A conditional random field model with additional codec used for
+-- data encoding.
+data CRF a b = CRF {
+    -- | The codec is used to transform data into internal representation,
+    -- where each observation and each label is represented by a unique
+    -- integer number.
+    codec :: Codec a b,
+    -- | The actual model, which is a map from 'Feature's to potentials.
+    model :: Model }
+
+instance (Ord a, Ord b, Binary a, Binary b) => Binary (CRF a b) where
+    put CRF{..} = put codec >> put model
+    get = CRF <$> get <*> get
+
+-- | Train the CRF using the stochastic gradient descent method.
+-- The resulting model will contain features extracted with
+-- the user supplied extraction function.
+-- You can use the functions provided by the "Data.CRF.Chain1.Feature.Present"
+-- and "Data.CRF.Chain1.Feature.Hidden" modules for this purpose.
+-- When the evaluation data 'IO' action is 'Just', the iterative
+-- training process will notify the user about the current accuracy
+-- on the evaluation part every full iteration over the training part.
+-- TODO: Accept custom r0 construction function.
+train
+    :: (Ord a, Ord b)
+    => SGD.SgdArgs                  -- ^ Args for SGD
+    -> IO [SentL a b]               -- ^ Training data 'IO' action
+    -> Maybe (IO [SentL a b])       -- ^ Maybe evalation data
+    -> (AVec Lb -> [(Xs, Ys)] -> [Feature])     -- ^ Feature selection
+    -> IO (CRF a b)                 -- ^ Resulting model
+train sgdArgs trainIO evalIO'Maybe extractFeats = do
+    hSetBuffering stdout NoBuffering
+    (_codec, trainData) <- mkCodec <$> trainIO
+    _r0 <- encodeLabels _codec . S.toList . unkSet <$> trainIO
+    evalDataM <- case evalIO'Maybe of
+        Just evalIO -> Just . encodeDataL _codec <$> evalIO
+        Nothing     -> return Nothing
+    let feats = extractFeats _r0 trainData
+        crf = (mkModel feats) { r0 = _r0 }
+    para <- SGD.sgdM sgdArgs
+        (notify sgdArgs crf trainData evalDataM)
+        (gradOn crf) (V.fromList trainData) (values crf)
+    return $ CRF _codec (crf { values = para })
+
+-- | Collect labels assigned to unknown words (with empty list
+-- of potential interpretations).
+unkSet :: Ord b => [SentL a b] -> S.Set b
+unkSet =
+    S.fromList . concatMap onSent
+  where
+    onSent = concatMap onWord
+    onWord word
+        | unknown (fst word)    = M.keys . unDist . snd $ word
+        | otherwise             = []
+
+gradOn :: Model -> SGD.Para -> (Xs, Ys) -> SGD.Grad
+gradOn crf para (xs, ys) = SGD.fromLogList $
+    [ (featToJustInt curr feat, L.fromPos val)
+    | (feat, val) <- featuresIn xs ys ] ++
+    [ (ix, L.fromNeg val)
+    | (FeatIx ix, val) <- expectedFeaturesIn curr xs ]
+  where
+    curr = crf { values = para }
+
+notify
+    :: SGD.SgdArgs -> Model -> [(Xs, Ys)] -> Maybe [(Xs, Ys)]
+    -> SGD.Para -> Int -> IO ()
+notify SGD.SgdArgs{..} crf trainData evalDataM para k 
+    | doneTotal k == doneTotal (k - 1) = putStr "."
+    | Just dataSet <- evalDataM = do
+        let x = accuracy (crf { values = para }) dataSet
+        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 = length trainData
diff --git a/Data/CRF/Chain1/Constrained/Util.hs b/Data/CRF/Chain1/Constrained/Util.hs
new file mode 100644
--- /dev/null
+++ b/Data/CRF/Chain1/Constrained/Util.hs
@@ -0,0 +1,12 @@
+module Data.CRF.Chain1.Constrained.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)
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 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-chain1-constrained.cabal b/crf-chain1-constrained.cabal
new file mode 100644
--- /dev/null
+++ b/crf-chain1-constrained.cabal
@@ -0,0 +1,71 @@
+name:               crf-chain1-constrained
+version:            0.1.0
+synopsis:           First-order, constrained, linear-chain conditional random fields
+description:
+    The library provides efficient implementation of the first-order,
+    linear-chain conditional random fields (CRFs) with position-wise
+    constraints imposed over label values.
+    .
+    It is strongly related to the simpler
+    <http://hackage.haskell.org/package/crf-chain1>
+    library where constraints are not taken into account and all
+    features which are not included in the CRF model are considered to have
+    probability of 0.  Here, on the other hand, such features do not influence
+    the overall probability of the (sentence, labels) pair - they are
+    assigned the default potential of 0.
+    .
+    Efficient algorithm for determining marginal probabilities of individual
+    labels is provided.
+    The tagging is performed with respect to marginal probabilities.
+--     The argmax algorithm (finding
+--     the most probable label sequence satisfying the given constraints)
+--     is less efficient, since we cannot use the sparse
+--     forward-backward recursions optimization.
+license:            BSD3
+license-file:       LICENSE
+cabal-version:      >= 1.6
+copyright:          Copyright (c) 2012 IPI PAN
+author:             Jakub Waszczuk
+maintainer:         waszczuk.kuba@gmail.com
+stability:          experimental
+category:           Math
+homepage:           https://github.com/kawu/crf-chain1-constrained
+build-type:         Simple
+
+library
+    build-depends:
+        base >= 4 && < 5
+      , containers
+      , vector
+      , array
+      , random
+      , parallel
+      , logfloat
+      , monad-codec >= 0.2 && < 0.3
+      , binary
+      , vector-binary
+      , data-lens
+      , sgd >= 0.2.1 && < 0.3
+
+    exposed-modules:
+        Data.CRF.Chain1.Constrained
+      , Data.CRF.Chain1.Constrained.Dataset.Internal
+      , Data.CRF.Chain1.Constrained.Dataset.External
+      , Data.CRF.Chain1.Constrained.Dataset.Codec
+      , Data.CRF.Chain1.Constrained.Feature
+      , Data.CRF.Chain1.Constrained.Feature.Present
+      , Data.CRF.Chain1.Constrained.Feature.Hidden
+      , Data.CRF.Chain1.Constrained.Model
+      , Data.CRF.Chain1.Constrained.Inference
+      , Data.CRF.Chain1.Constrained.Train
+
+    other-modules:
+        Data.CRF.Chain1.Constrained.DP
+      , Data.CRF.Chain1.Constrained.Util
+      , Data.CRF.Chain1.Constrained.Intersect
+        
+    ghc-options: -Wall -O2
+
+source-repository head
+    type: git
+    location: git://github.com/kawu/crf-chain1-constrained.git