diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2012, Grzegorz Chrupala
+
+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.
+
+    * Neither the name of Grzegorz Chrupala nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+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/NLP/SwiftLDA.hs b/NLP/SwiftLDA.hs
new file mode 100644
--- /dev/null
+++ b/NLP/SwiftLDA.hs
@@ -0,0 +1,367 @@
+{-# LANGUAGE NoMonomorphismRestriction 
+, BangPatterns 
+, DeriveGeneric
+ #-}
+-- | Latent Dirichlet Allocation
+--
+-- Imperative implementation of a collapsed Gibbs sampler for LDA. This
+-- library uses the topic modeling terminology (documents, words,
+-- topics), even though it is generic. For example if used for word
+-- class induction, replace documents with word types, words with
+-- features and topics with word classes.
+module NLP.SwiftLDA
+       ( -- * Samplers
+         pass
+       , passOne
+         -- * Datatypes
+       , LDA
+       , Doc
+       , D
+       , W
+       , Z
+       , Table2D
+       , Table1D
+         -- * Access model information
+       , Finalized (..)         
+         -- * Initialization and finalization
+       , initial
+       , finalize
+         -- * Querying evolving model
+       , docTopicWeights_
+       , priorDocTopicWeights_
+         -- * Querying finalized model
+       , docTopicWeights
+       , wordTopicWeights
+       )
+where       
+-- Standard libraries  
+import Prelude hiding (read, exponent)
+import Data.Array.ST
+import Data.STRef
+import Control.Applicative 
+import System.Random.MWC
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.ST
+import Data.Word
+import qualified Data.Vector.Unboxed         as U
+import qualified Data.Vector                 as V
+import qualified Data.IntMap                 as IntMap
+import qualified Data.List                   as List
+import qualified Data.Foldable               as Fold
+import GHC.Generics (Generic)
+import Debug.Trace
+
+-- Package modules
+import NLP.SwiftLDA.UnboxedMaybeVector ()
+
+type Array2D s = STUArray s (Int,Int) Double
+type Array1D s = STUArray s Int Double
+
+type Table2D = IntMap.IntMap Table1D
+type Table1D = IntMap.IntMap Double
+
+type D = Int
+type Z = Int  
+type W = Int
+type Doc = (D, U.Vector (W, Maybe Z))
+
+-- | Abstract type holding the settings and the state of the sampler
+data LDA s = 
+  LDA 
+  { _docTopics   :: !(STRef s (Array2D s))    -- ^ Document-topic counts
+  , _wordTopics  :: !(STRef s (Array2D s))    -- ^ Word-topic counts
+  , _topics      :: !(Array1D s)              -- ^ Topic counts
+  , _alphasum    :: !Double                   -- ^ alpha * K Dirichlet
+                                             -- parameter (topic
+                                             -- sparseness)
+  , _beta        :: !Double                   -- ^ beta Dirichlet
+                                             -- parameter (word
+                                             -- sparseness)
+  , _topicNum    :: !Int                      -- ^ Number of topics K
+  , _wSize       :: !(STRef s Int)            -- ^ Number of unique words
+  , weights     :: !(Array1D s)              -- ^ Current token weights
+  , weightSum   :: !(STRef s Double)         -- ^ Sum of current token weights
+  , gen         :: !(Gen (PrimState (ST s))) -- ^ Random generator
+  , _exponent   :: !(Maybe Double)
+  } 
+  
+data Finalized = 
+  Finalized 
+  { docTopics  :: !Table2D  -- ^ Document topic counts
+  , wordTopics :: !Table2D  -- ^ Word topic counts
+  , topics     :: !Table1D  -- ^ Topics counts
+  , topicDocs  :: !Table2D  -- ^ Inverse document-topic counts
+  , topicWords :: !Table2D  -- ^ Inverse word-topic counts
+  , alphasum   :: !Double   -- ^ alpha * K Dirichlet parameter (topic
+                                             -- sparseness)
+  , beta       :: !Double   -- ^ beta Dirichlet parameter (word
+                            -- sparseness)
+  , topicNum   :: !Int      -- ^ Number of topics K
+  , wSize      :: !Int      -- ^ Number of unique words
+  , exponent   :: !(Maybe Double)   -- ^ Learning rate exponent 
+  }  deriving (Generic)
+  
+-- | Create transparent immutable object holding model information
+-- from opaque internal representation
+finalize :: LDA s -> ST s Finalized
+finalize m = do
+  dt  <- read . _docTopics $ m
+  wt  <- read . _wordTopics $ m
+  dtf <- freezeArray2D dt
+  wtf <- freezeArray2D wt
+  tf  <- freezeArray1D (_topics m)
+  ws  <- read . _wSize $ m
+  return $! Finalized {
+      docTopics  = dtf
+    , wordTopics = wtf
+    , topics     = tf 
+    , topicDocs  = invert dtf
+    , topicWords = invert wtf
+    , alphasum   = _alphasum m
+    , beta       = _beta m
+    , topicNum   = _topicNum m 
+    , wSize      = ws
+    , exponent   = _exponent m
+    }
+  
+-- | Initial document index upper bound
+iDSIZE :: Int
+iDSIZE = 1000
+
+-- | Initial word index upper bound
+iWSIZE :: Int
+iWSIZE = 1000
+
+-- | @initial s k a b@ initializes model with @k@ topics, @a/k@ alpha
+-- hyperparameter, @b@ beta hyperparameter and random seed @s@
+initial :: U.Vector Word32 -> Int -> Double -> Double -> Maybe Double 
+           -> ST s (LDA s) 
+initial s k a b e = do           
+  dta <- newArray ((0,0),(iDSIZE, k-1)) 0
+  wta <- newArray ((0,0),(iWSIZE, k-1)) 0
+  LDA <$> 
+    new dta <*> 
+    new wta <*> 
+    newArray (0,k-1) 0 <*> 
+    pure a <*>
+    pure b <*>
+    pure k <*>
+    new 0 <*>
+    newArray (0,k-1) 0 <*>
+    new 0 <*>
+    initialize s <*>
+    pure e
+                   
+rho :: Double -> Int -> Double
+rho e t = 1 - (1 + fromIntegral t)**(-e)
+{-# INLINE rho #-}
+
+-- | @pass batch@ runs one pass of Gibbs sampling on documents in @batch@  
+pass :: Int -> LDA s -> V.Vector Doc -> ST s (V.Vector Doc)
+pass t m = V.mapM (passOne t m) 
+
+-- | Run a pass on a single doc  
+passOne :: Int -> LDA s -> Doc -> ST s Doc
+passOne t m doc@(!d, wz) = do
+  grow m doc
+  zs <- U.mapM one wz
+  return (d, U.zip (U.map fst wz) (U.map Just zs))
+  where r = maybe 1 (flip rho t) . _exponent $ m
+        one (w, mz) = do
+          case mz of
+            Nothing -> return ()
+            Just z   -> update (negate r) m d w z
+          !z <- randomZ m d w
+          update r m d w z
+          return z
+
+-- | Sample a random topic for doc d and word w        
+randomZ :: LDA s -> Int -> Int -> ST s Int
+randomZ m !d !w = do
+  wordTopicWeights_ m d w
+  !s <- read (weightSum m)
+  sample (weights m) s (gen m)
+
+  
+-- | @topicWeights m d w@ sets @weights m@ to the unnormalized probabilities of
+-- topics for word @w@ in document @d@ given LDA model @m@.  
+-- Each call overwrites current weights
+wordTopicWeights_ :: LDA s -> Int -> Int -> ST s ()
+wordTopicWeights_ m !d !w = do
+  let k  = _topicNum m
+      a  = _alphasum m / fromIntegral k
+      b  = _beta m
+  v  <- fromIntegral  <$> read (_wSize m)
+  dt <- read (_docTopics m)
+  wt <- read (_wordTopics m)
+  let ws = weights m
+  write (weightSum m) 0
+  (l,u) <- getBounds ws
+  let go !z | z > u = return ()
+      go !z = do
+        nzd <- readArray dt (d,z) 
+        nzw <- readArray wt (w,z)
+        nz  <- readArray (_topics m) z
+        let !n = (nzd + a) * (nzw + b) / (nz + v * b)
+        !s <- read (weightSum m)    
+        write (weightSum m) (s+n)
+        writeArray ws z n
+        go (z+1)
+  go l      
+  
+-- Weight calc on Finalized
+-- | @docTopicWeights m doc@ returns unnormalized topic probabilities
+-- for document doc given LDA model @m@
+docTopicWeights :: Finalized -> Doc -> U.Vector Double
+docTopicWeights m (d, ws) = 
+    U.accumulate (+) (U.replicate (topicNum m) 0)
+  . U.concatMap (U.indexed . wordTopicWeights m d)
+  . U.map fst 
+  $ ws
+{-# INLINE docTopicWeights #-}
+    
+priorDocTopicWeights_ :: LDA s -> D -> ST s (U.Vector Double)
+priorDocTopicWeights_ m d = do
+  grow m (d, U.empty)
+  dt <- read (_docTopics m)
+  ((_,0),(_,u)) <- getBounds dt
+  U.generateM (u+1) (\z -> readArray dt (d,z))
+  
+docTopicWeights_ :: LDA s -> Doc -> ST s (U.Vector Double)
+docTopicWeights_ m doc@(d, ws) = do
+  grow m doc
+  (0,u) <- getBounds (weights m)
+  let r = U.replicate (_topicNum m) 0
+  let one w = do
+        wordTopicWeights_ m d w
+        U.generateM (u+1) (readArray (weights m))
+  U.foldM' (\z w -> do y <- one w 
+                       return $! U.zipWith (+) z y) r 
+    . U.map fst 
+    $ ws
+-- | @topicWeights m d w@ returns the unnormalized probabilities of
+-- topics for word @w@ in document @d@ given LDA model @m@.
+wordTopicWeights :: Finalized -> D -> W -> U.Vector Double
+wordTopicWeights m d w =
+  let k = topicNum m
+      a = alphasum m / fromIntegral k
+      b = beta m
+      dt = IntMap.findWithDefault IntMap.empty d . docTopics $ m
+      wt = IntMap.findWithDefault IntMap.empty w . wordTopics $ m
+      v = fromIntegral . wSize $ m
+      weights = [   (count z dt + a)     -- n(z,d) + alpha 
+                  * (count z wt + b)     -- n(z,w) + beta
+                  * (1/(count z (topics m) + v * b)) 
+                      -- n(.,w) + V * beta
+                  | z <- [0..k-1] ]
+  in U.fromList weights
+{-# INLINE wordTopicWeights #-}
+
+-- | Update counts in the model corresponding to given doc, word and topic
+update :: Double -> LDA s -> Int -> Int -> Int -> ST s ()  
+update c m d w z = do
+  dt  <- read (_docTopics m)
+  wt  <- read (_wordTopics m)
+  wsz <- read (_wSize m) ; write (_wSize m) (max (w+1) wsz)
+  nz  <- readArray (_topics m) z ; writeArray (_topics m) z (nz+c)
+  nzd <- readArray dt (d,z)     ; writeArray dt (d,z) (nzd+c)
+  nzw <- readArray wt (w,z)     ; writeArray wt (w,z) (nzw+c)
+  
+-- | @grow m doc@ grows the @docTopic m@ and @wordTopic m@ tables as necessary
+-- according to content of @doc@
+grow :: LDA s -> Doc -> ST s ()
+grow m (d, wz) = do
+  let w = if U.null wz then 0 else U.maximum  (U.map fst wz)
+  dt <- read (_docTopics m)  ; (_,(d_max,_)) <- getBounds dt 
+  wt <- read (_wordTopics m) ; (_,(w_max,_)) <- getBounds wt
+  when (d > d_max) (do dt' <- resize dt
+                       write (_docTopics m)  dt')
+  when (w > w_max) (do wt' <- resize wt
+                       write (_wordTopics m) wt')
+  
+-- | @resize table@ creates a new array with double the size of the
+-- first component of the upper bound of @arr@ and copies to content
+-- of @arr@ to it.
+resize :: Array2D s -> ST s (Array2D s)
+resize a = do
+  bs@((l1,l2),(u1_old,u2)) <- getBounds a
+  trace (show bs) () `seq` return ()
+  let u1 = u1_old * 2
+      bs' = ((l1,l2),(u1,u2))
+  b <- newArray bs' 0
+  let copy !i = do
+        v <- readArray a i
+        writeArray b i v
+  mapM_ copy (range bs)         
+  return b
+
+-- | @sample ws s g@ draws a random topic according to weights @ws@
+-- which should sum up to @s@
+sample :: Array1D s -> Double -> Gen s -> ST s Int  
+sample ws s g = do
+  !r <- uniformR (0,s) g
+  findEvent r ws
+  
+-- | @findEvent r ws@ converts weights to cumulative weights, and
+-- finds the index whose cumulative weight is >= r.
+findEvent :: Double -> Array1D s -> ST s Int  
+findEvent !r ws = do
+  (l,u) <- getBounds ws
+  let go !i !_n | i > u = return (i-1) 
+      go !i !n  | n > 0.0 = do v <- readArray ws i
+                               go (i+1) (n-v)
+                | otherwise = return (i-1)  
+  go l r
+
+
+read :: STRef s a -> ST s a             
+read = readSTRef
+{-# INLINE read #-}
+
+write :: STRef s a -> a -> ST s ()
+write = writeSTRef
+{-# INLINE write #-}
+
+new :: a -> ST s (STRef s a) 
+new = newSTRef  
+{-# INLINE new #-}
+
+-- | Swap the order of keys on Table2D
+invert :: Table2D -> Table2D
+invert outer = 
+  List.foldl' (\z (k,k',v) -> upd v z k k')  IntMap.empty 
+  [ (k',k,v)
+    | (k, inner) <- IntMap.toList outer
+    , (k', v) <- IntMap.toList inner ]
+{-# INLINE invert #-}  
+
+-- | Increment table m by c at key (k,k')
+upd :: Double -> Table2D -> Int -> Int -> Table2D
+upd c m k k' = IntMap.insertWith' (flip (IntMap.unionWith (+)))
+                               k 
+                               (IntMap.singleton k' c)
+                               m
+
+{-# INLINE upd #-}
+
+freezeArray2D :: Array2D s -> ST s Table2D
+freezeArray2D a = do
+  bs <- getBounds a
+  Fold.foldlM f IntMap.empty (range bs)
+  where f z ind@(!i,!i') = do 
+          !v <- readArray a ind
+          if v > 0 
+            then return $! upd v z i i'
+            else return $! z
+                 
+freezeArray1D :: Array1D s -> ST s Table1D
+freezeArray1D a = IntMap.fromList . filter ((>0) . snd) <$> getAssocs a
+
+
+count :: Int -> IntMap.IntMap Double -> Double
+count z t = case IntMap.findWithDefault 0 z t of
+        n | n < 0 -> error "NLP.SwiftLDA.count: negative count"
+        n -> n
+{-# INLINE count #-}     
+        
diff --git a/NLP/SwiftLDA/UnboxedMaybeVector.hs b/NLP/SwiftLDA/UnboxedMaybeVector.hs
new file mode 100644
--- /dev/null
+++ b/NLP/SwiftLDA/UnboxedMaybeVector.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE Rank2Types, MultiParamTypeClasses, FlexibleContexts,
+             TypeFamilies, ScopedTypeVariables, BangPatterns,  
+             TupleSections  #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+module NLP.SwiftLDA.UnboxedMaybeVector ()
+where
+import qualified Data.Vector.Generic         as G
+import qualified Data.Vector.Generic.Mutable as M
+import Data.Vector.Unboxed.Base (MVector, Vector, Unbox)
+
+import Control.Monad ( liftM )
+
+newtype instance MVector s (Maybe a) = MV_Maybe (MVector s (a,Bool))
+newtype instance Vector    (Maybe a) = V_Maybe  (Vector    (a,Bool))
+
+instance (Num a, Unbox a) => Unbox (Maybe a)
+
+nothing :: (Num a, Unbox a) => a
+nothing = 0
+{-# INLINE nothing #-}
+
+instance (Num a, Unbox a) => M.MVector MVector (Maybe a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  basicLength (MV_Maybe v) = M.basicLength v
+  basicUnsafeSlice i n (MV_Maybe v) = MV_Maybe $ M.basicUnsafeSlice i n v
+  basicOverlaps (MV_Maybe v1) (MV_Maybe v2) = M.basicOverlaps v1 v2
+  basicUnsafeNew n = MV_Maybe `liftM` M.basicUnsafeNew n
+  basicUnsafeRead (MV_Maybe v) i = fromPair `liftM` M.basicUnsafeRead v i
+  basicUnsafeWrite (MV_Maybe v) i mx = M.basicUnsafeWrite v i (maybe (nothing,False) (,True) mx)
+  
+instance (Num a, Unbox a) => G.Vector Vector (Maybe a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  basicLength (V_Maybe v) = G.basicLength v
+  basicUnsafeFreeze (MV_Maybe v) = V_Maybe `liftM` G.basicUnsafeFreeze v
+  basicUnsafeThaw (V_Maybe v) = MV_Maybe `liftM` G.basicUnsafeThaw v
+  basicUnsafeSlice i n (V_Maybe v) = V_Maybe $ G.basicUnsafeSlice i n v
+  basicUnsafeIndexM (V_Maybe v) i
+                = fromPair `liftM` G.basicUnsafeIndexM v i
+                  
+fromPair :: (Unbox a) => (a, Bool) -> Maybe a                  
+fromPair (_, False) = Nothing
+fromPair (x, True)  = Just x
+{-# INLINE fromPair #-}
diff --git a/README b/README
new file mode 100644
--- /dev/null
+++ b/README
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/swift-lda.cabal b/swift-lda.cabal
new file mode 100644
--- /dev/null
+++ b/swift-lda.cabal
@@ -0,0 +1,33 @@
+Name:                swift-lda
+Version:             0.4.0
+Synopsis:            Online sampler for Latent Dirichlet Allocation
+Description:     Online Gibbs sampler for Latent Dirichlet Allocation. 
+                 LDA is a generative admixture model frequently used 
+                 for topic modeling and other applications. The primary
+                 goal of this implementation is to be used for probabilistic 
+                 soft word class induction.
+                 The sampler can be used in an online as well as batch mode.
+                 This package uses an imperative implementation in the ST monad.
+Homepage:            https://bitbucket.org/gchrupala/colada
+License:             BSD3
+License-file:        LICENSE
+Author:              Grzegorz Chrupała <pitekus@gmail.com>
+Maintainer:          Grzegorz Chrupała <pitekus@gmail.com>
+Category:            Natural Language Processing
+Extra-source-files:  README
+Build-type:          Simple
+Cabal-version:       >=1.2
+
+
+Library
+   Exposed-modules: NLP.SwiftLDA
+                  , NLP.SwiftLDA.UnboxedMaybeVector
+   Build-depends: base >= 3 && < 5
+                , vector >= 0.9
+                , primitive >= 0.4
+                , mwc-random >= 0.12
+                , array >= 0.3
+                , ghc-prim >= 0.2
+                , containers >= 0.4
+  Other-modules: 
+  GHC-options: -O2