diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,5 @@
+
+Version 0.1
+-----
+
+* Initial version.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2016 David Sorokin <david.sorokin@gmail.com>
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. 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.
+
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/Simulation/Aivika/Lattice.hs b/Simulation/Aivika/Lattice.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice.hs
@@ -0,0 +1,27 @@
+
+-- |
+-- Module     : Simulation.Aivika.Lattice
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- This module re-exports the library functionality related to nested computations
+-- that can be run within lattice nodes.
+--
+module Simulation.Aivika.Lattice
+       (-- * Modules
+        module Simulation.Aivika.Lattice.LIO,
+        module Simulation.Aivika.Lattice.Estimate,
+        module Simulation.Aivika.Lattice.Event,
+        module Simulation.Aivika.Lattice.Generator,
+        module Simulation.Aivika.Lattice.QueueStrategy,
+        module Simulation.Aivika.Lattice.Ref.Base) where
+
+import Simulation.Aivika.Lattice.LIO
+import Simulation.Aivika.Lattice.Estimate
+import Simulation.Aivika.Lattice.Event
+import Simulation.Aivika.Lattice.Generator
+import Simulation.Aivika.Lattice.QueueStrategy
+import Simulation.Aivika.Lattice.Ref.Base
diff --git a/Simulation/Aivika/Lattice/Estimate.hs b/Simulation/Aivika/Lattice/Estimate.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Estimate.hs
@@ -0,0 +1,147 @@
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Estimate
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 8.0.1
+--
+-- The module defines the 'Estimate' monad transformer which is destined for estimating
+-- computations within lattice nodes. Such computations are separated from the 'Event'
+-- computations. An idea is that the forward-traversing 'Event' computations provide with
+-- something that can be observed, while the backward-traversing 'Estimate' computations
+-- estimate the received information.
+--
+module Simulation.Aivika.Lattice.Estimate
+       (-- * Estimate Monad
+        Estimate,
+        EstimateLift(..),
+        runEstimateInStartTime,
+        estimateTime,
+        latticeTimeStep,
+        -- * Computations within Lattice
+        foldEstimate,
+        memoEstimate,
+        estimateUpSide,
+        estimateDownSide,
+        estimateFuture,
+        -- * Error Handling
+        catchEstimate,
+        finallyEstimate,
+        throwEstimate,
+        -- * Debugging
+        traceEstimate) where
+
+import Simulation.Aivika.Trans.Internal.Types
+import Simulation.Aivika.Trans.Simulation
+import Simulation.Aivika.Trans.Parameter
+import Simulation.Aivika.Trans.Ref
+import Simulation.Aivika.Trans.Observable
+import Simulation.Aivika.Lattice.Internal.Estimate
+import Simulation.Aivika.Lattice.Internal.LIO
+import qualified Simulation.Aivika.Lattice.Internal.Ref as R
+
+-- | Estimate the computation in the lattice nodes.
+memoEstimate :: (Estimate LIO a -> Estimate LIO a)
+                -- ^ estimate in the intermediate time point of the lattice
+                -> Estimate LIO a
+                -- ^ estimate in the final time point of the lattice or beyond it
+                -> Simulation LIO (Estimate LIO a)
+memoEstimate f m =
+  do r  <- R.newRef Nothing
+     t2 <- liftParameter stoptime
+     let loop =
+           Estimate $ \p ->
+           do b <- R.readRef0 r
+              case b of
+                Just a  -> return a
+                Nothing ->
+                  if pointTime p >= t2
+                  then do a <- invokeEstimate p m
+                          R.writeRef0 r (Just a)
+                          return a
+                  else do a <- invokeEstimate p (f loop)
+                          R.writeRef0 r (Just a)
+                          return a
+     return loop
+
+-- | Estimate the computation in the up side node of the lattice,
+-- where 'latticeTimeIndex' is increased by 1 but 'latticeMemberIndex' remains the same.
+--
+-- It is merely equivalent to the following definition:
+--
+-- @estimateUpSide = estimateFuture 1 0@
+--
+estimateUpSide :: Estimate LIO a -> Estimate LIO a
+estimateUpSide m =
+  Estimate $ \p ->
+  LIO $ \ps ->
+  do let ps' = upSideLIOParams ps
+         r   = pointRun p
+     p' <- invokeLIO ps' $
+           invokeParameter r
+           latticePoint
+     invokeLIO ps' $
+       invokeEstimate p' m
+
+-- | Estimate the computation in the down side node of the lattice,
+-- where the both 'latticeTimeIndex' and 'latticeMemberIndex' are increased by 1.
+--
+-- It is merely equivalent to the following definition:
+--
+-- @estimateDownSide = estimateFuture 1 1@
+--
+estimateDownSide :: Estimate LIO a -> Estimate LIO a
+estimateDownSide m =
+  Estimate $ \p ->
+  LIO $ \ps ->
+  do let ps' = downSideLIOParams ps
+         r   = pointRun p
+     p' <- invokeLIO ps' $
+           invokeParameter r
+           latticePoint
+     invokeLIO ps' $
+       invokeEstimate p' m
+
+-- | Estimate the computation in the shifted lattice node, where the first parameter
+-- specifies the positive 'latticeTimeIndex' shift, but the second parameter
+-- specifies the 'latticeMemberIndex' shift af any sign.
+--
+-- It allows looking into the future computations. The lattice is constructed in such a way
+-- that we can define the past 'Estimate' computation in terms of the future @Estimate@
+-- computation. That is the point.
+--
+-- Regarding the 'Event' computation, a quite opposite rule is true. The future @Event@ computation
+-- depends on the past @Event@ computation. But we can update 'Ref' references within
+-- the corresponding discrete event simulation and then read them within the @Estimate@
+-- computation, because @Ref@ is 'Observable'.
+estimateFuture :: Int
+                  -- ^ a positive shift of the lattice time index
+                  -> Int
+                  -- ^ a shift of the lattice member index
+                  -> Estimate LIO a
+                  -- ^ the source computation
+                  -> Estimate LIO a
+estimateFuture di dk m =
+  Estimate $ \p ->
+  LIO $ \ps ->
+  do let ps' = shiftLIOParams di dk ps
+         r   = pointRun p
+     p' <- invokeLIO ps' $
+           invokeParameter r
+           latticePoint
+     invokeLIO ps' $
+       invokeEstimate p' m
+
+-- | Fold the estimation of the specified computation.
+foldEstimate :: (a -> a -> Estimate LIO a)
+                -- ^ reduce in the intermediate nodes of the lattice
+                -> Estimate LIO a
+                -- ^ estimate the computation in the final time point and beyond it
+                -> Simulation LIO (Estimate LIO a)
+foldEstimate f = memoEstimate g
+  where g m =
+          do a1 <- estimateUpSide m
+             a2 <- estimateDownSide m
+             f a1 a2
diff --git a/Simulation/Aivika/Lattice/Event.hs b/Simulation/Aivika/Lattice/Event.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Event.hs
@@ -0,0 +1,18 @@
+
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Event
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- The module defines an event queue, where 'LIO' is an instance of 'EventQueueing'.
+--
+module Simulation.Aivika.Lattice.Event () where
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice.Internal.Event
+import Simulation.Aivika.Lattice.Internal.LIO
diff --git a/Simulation/Aivika/Lattice/Generator.hs b/Simulation/Aivika/Lattice/Generator.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Generator.hs
@@ -0,0 +1,125 @@
+
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Generator
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- Here is defined a random number generator,
+-- where 'LIO' is an instance of 'MonadGenerator'.
+--
+module Simulation.Aivika.Lattice.Generator () where
+
+import Control.Monad
+import Control.Monad.Trans
+
+import System.Random
+
+import Data.IORef
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Trans.Generator.Primitive
+import Simulation.Aivika.Lattice.Internal.LIO
+
+instance MonadGenerator LIO where
+
+  data Generator LIO =
+    GeneratorLIO { generator01 :: LIO Double,
+                   -- ^ the generator of uniform numbers from 0 to 1
+                   generatorNormal01 :: LIO Double
+                   -- ^ the generator of normal numbers with mean 0 and variance 1
+                 }
+
+  generateUniform = generateUniform01 . generator01
+
+  generateUniformInt = generateUniformInt01 . generator01
+
+  generateTriangular = generateTriangular01 . generator01
+
+  generateNormal = generateNormal01 . generatorNormal01
+
+  generateLogNormal = generateLogNormal01 . generatorNormal01
+
+  generateExponential = generateExponential01 . generator01
+
+  generateErlang = generateErlang01 . generator01
+
+  generatePoisson = generatePoisson01 . generator01
+
+  generateBinomial = generateBinomial01 . generator01
+
+  generateGamma g = generateGamma01 (generatorNormal01 g) (generator01 g)
+
+  generateBeta g = generateBeta01 (generatorNormal01 g) (generator01 g)
+
+  generateWeibull = generateWeibull01 . generator01
+
+  generateDiscrete = generateDiscrete01 . generator01
+
+  newGenerator tp =
+    case tp of
+      SimpleGenerator ->
+        liftIO newStdGen >>= newRandomGenerator
+      SimpleGeneratorWithSeed x ->
+        error "Unsupported generator type SimpleGeneratorWithSeed: newGenerator"
+      CustomGenerator g ->
+        g
+      CustomGenerator01 g ->
+        newRandomGenerator01 g
+
+  newRandomGenerator g = 
+    do r <- liftIO $ newIORef g
+       let g01 = do g <- liftIO $ readIORef r
+                    let (x, g') = random g
+                    liftIO $ writeIORef r g'
+                    return x
+       newRandomGenerator01 g01
+
+  newRandomGenerator01 g01 =
+    do gNormal01 <- newNormalGenerator01 g01
+       return GeneratorLIO { generator01 = g01,
+                             generatorNormal01 = gNormal01 }
+
+-- | Create a normal random number generator with mean 0 and variance 1
+-- by the specified generator of uniform random numbers from 0 to 1.
+newNormalGenerator01 :: LIO Double
+                        -- ^ the generator
+                        -> LIO (LIO Double)
+newNormalGenerator01 g =
+  do nextRef <- liftIO $ newIORef 0.0
+     flagRef <- liftIO $ newIORef False
+     xi1Ref  <- liftIO $ newIORef 0.0
+     xi2Ref  <- liftIO $ newIORef 0.0
+     psiRef  <- liftIO $ newIORef 0.0
+     let loop =
+           do psi <- liftIO $ readIORef psiRef
+              if (psi >= 1.0) || (psi == 0.0)
+                then do g1 <- g
+                        g2 <- g
+                        let xi1 = 2.0 * g1 - 1.0
+                            xi2 = 2.0 * g2 - 1.0
+                            psi = xi1 * xi1 + xi2 * xi2
+                        liftIO $ writeIORef xi1Ref xi1
+                        liftIO $ writeIORef xi2Ref xi2
+                        liftIO $ writeIORef psiRef psi
+                        loop
+                else liftIO $ writeIORef psiRef $ sqrt (- 2.0 * log psi / psi)
+     return $
+       do flag <- liftIO $ readIORef flagRef
+          if flag
+            then do liftIO $ writeIORef flagRef False
+                    liftIO $ readIORef nextRef
+            else do liftIO $ writeIORef xi1Ref 0.0
+                    liftIO $ writeIORef xi2Ref 0.0
+                    liftIO $ writeIORef psiRef 0.0
+                    loop
+                    xi1 <- liftIO $ readIORef xi1Ref
+                    xi2 <- liftIO $ readIORef xi2Ref
+                    psi <- liftIO $ readIORef psiRef
+                    liftIO $ writeIORef flagRef True
+                    liftIO $ writeIORef nextRef $ xi2 * psi
+                    return $ xi1 * psi
diff --git a/Simulation/Aivika/Lattice/Internal/Estimate.hs b/Simulation/Aivika/Lattice/Internal/Estimate.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Internal/Estimate.hs
@@ -0,0 +1,164 @@
+
+{-# LANGUAGE RecursiveDo, MultiParamTypeClasses, FlexibleInstances #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Internal.Estimate
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 8.0.1
+--
+-- The module defines the 'Estimate' monad transformer which is destined for estimating
+-- computations within lattice nodes. Such computations are separated from the 'Event'
+-- computations. An idea is that the forward-traversing 'Event' computations provide with
+-- something that can be observed, while the backward-traversing 'Estimate' computations
+-- estimate the received information.
+--
+module Simulation.Aivika.Lattice.Internal.Estimate
+       (-- * Estimate Monad
+        Estimate(..),
+        EstimateLift(..),
+        invokeEstimate,
+        runEstimateInStartTime,
+        estimateTime,
+        -- * Error Handling
+        catchEstimate,
+        finallyEstimate,
+        throwEstimate,
+        -- * Debugging
+        traceEstimate) where
+
+import Control.Exception
+import Control.Monad
+import Control.Monad.Trans
+import Control.Monad.Fix
+import Control.Applicative
+
+import Debug.Trace (trace)
+
+import Simulation.Aivika.Trans.Exception
+import Simulation.Aivika.Trans.Comp
+import Simulation.Aivika.Trans.DES
+import Simulation.Aivika.Trans.Internal.Types
+import Simulation.Aivika.Trans.Parameter
+import Simulation.Aivika.Trans.Dynamics
+import Simulation.Aivika.Trans.Event
+import Simulation.Aivika.Lattice.Internal.LIO
+
+-- | A value in the 'Estimate' monad transformer represents something
+-- that can be estimated within lattice nodes.
+newtype Estimate m a = Estimate (Point m -> m a)
+
+-- | Invoke the 'Estimate' computation.
+invokeEstimate :: Point m -> Estimate m a -> m a
+{-# INLINE invokeEstimate #-}
+invokeEstimate p (Estimate m) = m p
+
+instance Monad m => Monad (Estimate m) where
+
+  {-# INLINE return #-}
+  return a = Estimate $ \p -> return a
+
+  {-# INLINE (>>=) #-}
+  (Estimate m) >>= k =
+    Estimate $ \p -> 
+    do a <- m p
+       let Estimate m' = k a
+       m' p
+
+instance Functor m => Functor (Estimate m) where
+  
+  {-# INLINE fmap #-}
+  fmap f (Estimate x) = Estimate $ \p -> fmap f $ x p
+
+instance Applicative m => Applicative (Estimate m) where
+  
+  {-# INLINE pure #-}
+  pure = Estimate . const . pure
+  
+  {-# INLINE (<*>) #-}
+  (Estimate x) <*> (Estimate y) = Estimate $ \p -> x p <*> y p
+
+instance MonadTrans Estimate where
+
+  {-# INLINE lift #-}
+  lift = Estimate . const
+
+instance MonadIO m => MonadIO (Estimate m) where
+  
+  {-# INLINE liftIO #-}
+  liftIO = Estimate . const . liftIO
+
+instance MonadFix m => MonadFix (Estimate m) where
+
+  {-# INLINE mfix #-}
+  mfix f = 
+    Estimate $ \p ->
+    do { rec { a <- invokeEstimate p (f a) }; return a }
+
+instance Monad m => MonadCompTrans Estimate m where
+
+  {-# INLINE liftComp #-}
+  liftComp = Estimate . const
+
+-- | A type class to lift the 'Estimate' computations into other computations.
+class EstimateLift t m where
+  
+  -- | Lift the specified 'Estimate' computation into another computation.
+  liftEstimate :: Estimate m a -> t m a
+
+instance Monad m => EstimateLift Estimate m where
+  
+  {-# INLINE liftEstimate #-}
+  liftEstimate = id
+
+instance Monad m => ParameterLift Estimate m where
+
+  {-# INLINE liftParameter #-}
+  liftParameter (Parameter x) = Estimate $ x . pointRun
+
+-- | Exception handling within 'Estimate' computations.
+catchEstimate :: (MonadException m, Exception e) => Estimate m a -> (e -> Estimate m a) -> Estimate m a
+{-# INLINABLE catchEstimate #-}
+catchEstimate (Estimate m) h =
+  Estimate $ \p -> 
+  catchComp (m p) $ \e ->
+  let Estimate m' = h e in m' p
+                           
+-- | A computation with finalization part like the 'finally' function.
+finallyEstimate :: MonadException m => Estimate m a -> Estimate m b -> Estimate m a
+{-# INLINABLE finallyEstimate #-}
+finallyEstimate (Estimate m) (Estimate m') =
+  Estimate $ \p ->
+  finallyComp (m p) (m' p)
+
+-- | Like the standard 'throw' function.
+throwEstimate :: (MonadException m, Exception e) => e -> Estimate m a
+{-# INLINABLE throwEstimate #-}
+throwEstimate e =
+  Estimate $ \p ->
+  throwComp e
+
+-- | Run the 'Estimate' computation in the start time and return the estimate.
+runEstimateInStartTime :: MonadDES m => Estimate m a -> Simulation m a
+{-# INLINE runEstimateInStartTime #-}
+runEstimateInStartTime (Estimate m) = runEventInStartTime (Event m)
+
+-- | Like 'time' estimates the current modeling time.
+estimateTime :: MonadDES m => Estimate m Double
+{-# INLINE estimateTime #-}
+estimateTime = Estimate $ return . pointTime
+
+-- | Show the debug message with the current simulation time and lattice node indices.
+traceEstimate :: String -> Estimate LIO a -> Estimate LIO a
+{-# INLINABLE traceEstimate #-}
+traceEstimate message m =
+  Estimate $ \p ->
+  LIO $ \ps ->
+  trace ("t = " ++ show (pointTime p) ++
+         ", lattice time index = " ++ show (lioTimeIndex ps) ++
+         ", lattice member index = " ++ show (lioMemberIndex ps) ++
+         ": " ++ message) $
+  invokeLIO ps $
+  invokeEstimate p m
diff --git a/Simulation/Aivika/Lattice/Internal/Event.hs b/Simulation/Aivika/Lattice/Internal/Event.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Internal/Event.hs
@@ -0,0 +1,206 @@
+
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Internal.Event
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- The module defines an event queue, where 'LIO' is an instance of 'EventQueueing'.
+-- Also it defines basic functions for running nested computations within lattice nodes.
+--
+module Simulation.Aivika.Lattice.Internal.Event
+       (estimateRef) where
+
+import Data.IORef
+
+import Control.Monad
+import Control.Monad.Trans
+
+import qualified Simulation.Aivika.PriorityQueue.Pure as PQ
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Trans.Internal.Types
+
+import Simulation.Aivika.Lattice.Internal.LIO
+import Simulation.Aivika.Lattice.Internal.Estimate
+import qualified Simulation.Aivika.Lattice.Internal.Ref as R
+
+-- | An implementation of the 'EventQueueing' type class.
+instance EventQueueing LIO where
+
+  -- | The event queue type.
+  data EventQueue LIO =
+    EventQueueLIO { queuePQ :: R.Ref (PQ.PriorityQueue (Point LIO -> LIO ())),
+                    -- ^ the underlying priority queue
+                    queueBusy :: IORef Bool,
+                    -- ^ whether the queue is currently processing events
+                    queueTime :: R.Ref Double
+                    -- ^ the actual time of the event queue
+                  }
+
+  newEventQueue specs =
+    do f  <- liftIO $ newIORef False
+       t  <- R.newRef0 (spcStartTime specs)
+       pq <- R.newRef0 PQ.emptyQueue
+       return EventQueueLIO { queuePQ   = pq,
+                              queueBusy = f,
+                              queueTime = t }
+
+  enqueueEvent t (Event m) =
+    Event $ \p ->
+    let pq = queuePQ $ runEventQueue $ pointRun p
+    in invokeEvent p $
+       R.modifyRef pq $ \x -> PQ.enqueue x t m
+
+  runEventWith processing (Event e) =
+    Dynamics $ \p ->
+    LIO $ \ps ->
+    do invokeLIO ps $
+         invokeDynamics p $
+         processEvents processing
+       invokeLIO ps $
+         e p
+
+  eventQueueCount =
+    Event $ \p ->
+    let pq = queuePQ $ runEventQueue $ pointRun p
+    in invokeEvent p $
+       fmap PQ.queueCount $ R.readRef pq
+
+-- | Process the pending events.
+processPendingEventsCore :: Bool -> Dynamics LIO ()
+processPendingEventsCore includingCurrentEvents = Dynamics r where
+  r p =
+    LIO $ \ps ->
+    do let q = runEventQueue $ pointRun p
+           f = queueBusy q
+       f' <- readIORef f
+       if f'
+         then error $
+              "Detected an event loop, which may indicate to " ++
+              "a logical error in the model: processPendingEventsCore"
+         else do writeIORef f True
+                 invokeLIO ps $
+                   invokeDynamics p $
+                   processPendingEventsUnsafe includingCurrentEvents
+                 writeIORef f False
+
+-- | Process the pending events in unsafe manner.
+processPendingEventsUnsafe :: Bool -> Dynamics LIO ()
+processPendingEventsUnsafe includingCurrentEvents = Dynamics r where
+  r p =
+    LIO $ \ps ->
+    let q = runEventQueue $ pointRun p
+    in call q p ps
+  call q p ps =
+    do let pq = queuePQ q
+           r  = pointRun p
+       f <- invokeLIO ps $
+            fmap PQ.queueNull $ R.readRef0 pq
+       unless f $
+         do (t2, c2) <- invokeLIO ps $
+                        fmap PQ.queueFront $ R.readRef0 pq
+            let t = queueTime q
+            t' <- invokeLIO ps $
+                  R.readRef0 t
+            when (t2 < t') $ 
+              error "The time value is too small: processPendingEventsCore"
+            when ((t2 < pointTime p) ||
+                  (includingCurrentEvents && (t2 == pointTime p))) $
+              do let sc = pointSpecs p
+                     t0 = spcStartTime sc
+                     dt = spcDT sc
+                     n2 = fromIntegral $ floor ((t2 - t0) / dt)
+                     p2 = p { pointTime = t2,
+                              pointIteration = n2,
+                              pointPhase = -1 }
+                 invokeLIO ps $
+                   R.writeRef0 t t2
+                 invokeLIO ps $
+                   R.defineTopRef0_ pq
+                 invokeLIO ps $
+                   R.modifyRef0 pq PQ.dequeue
+                 invokeLIO ps $
+                   c2 p2
+                 call q p ps
+
+-- | Process the pending events synchronously, i.e. without past.
+processPendingEvents :: Bool -> Dynamics LIO ()
+processPendingEvents includingCurrentEvents = Dynamics r where
+  r p =
+    LIO $ \ps ->
+    do let q = runEventQueue $ pointRun p
+           t = queueTime q
+       t' <- invokeLIO ps $
+             invokeEvent p $
+             R.readRef t
+       if pointTime p < t'
+         then error $
+              "The current time is less than " ++
+              "the time in the queue: processPendingEvents"
+         else invokeLIO ps $
+              invokeDynamics p $
+              processPendingEventsCore includingCurrentEvents
+
+-- | A memoized value.
+processEventsIncludingCurrent :: Dynamics LIO ()
+processEventsIncludingCurrent = processPendingEvents True
+
+-- | A memoized value.
+processEventsIncludingEarlier :: Dynamics LIO ()
+processEventsIncludingEarlier = processPendingEvents False
+
+-- | A memoized value.
+processEventsIncludingCurrentCore :: Dynamics LIO ()
+processEventsIncludingCurrentCore = processPendingEventsCore True
+
+-- | A memoized value.
+processEventsIncludingEarlierCore :: Dynamics LIO ()
+processEventsIncludingEarlierCore = processPendingEventsCore True
+
+-- | Process the events.
+processEvents :: EventProcessing -> Dynamics LIO ()
+processEvents CurrentEvents = processEventsIncludingCurrent
+processEvents EarlierEvents = processEventsIncludingEarlier
+processEvents CurrentEventsOrFromPast = processEventsIncludingCurrentCore
+processEvents EarlierEventsOrFromPast = processEventsIncludingEarlierCore
+
+-- | Initialize the event queue in the current lattice node if required.
+initEventQueue :: Event LIO ()
+initEventQueue =
+  Event $ \p ->
+  LIO $ \ps ->
+  do let pq = queuePQ $ runEventQueue r
+         r  = pointRun p
+     f <- invokeLIO ps $
+          R.topRefDefined0 pq
+     unless f $
+       do case parentLIOParams ps of
+            Nothing  -> error "The root must be initialized: initEventQueue"
+            Just ps' ->
+              do p' <- invokeLIO ps' $
+                       invokeParameter r
+                       latticePoint
+                 invokeLIO ps' $
+                   invokeEvent p'
+                   initEventQueue
+          invokeLIO ps $
+            R.defineTopRef0_ pq
+          invokeLIO ps $
+            invokeDynamics p $
+            processPendingEventsUnsafe True
+
+-- | Estimate the specified reference.
+estimateRef :: R.Ref a -> Estimate LIO a
+estimateRef r =
+  Estimate $ \p ->
+  LIO $ \ps ->
+  do invokeLIO ps $
+       invokeEvent p
+       initEventQueue
+     invokeLIO ps $
+       R.readRef0 r
diff --git a/Simulation/Aivika/Lattice/Internal/LIO.hs b/Simulation/Aivika/Lattice/Internal/LIO.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Internal/LIO.hs
@@ -0,0 +1,193 @@
+
+{-# LANGUAGE RecursiveDo #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Internal.LIO
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- This module defines the 'LIO' computation.
+--
+module Simulation.Aivika.Lattice.Internal.LIO
+       (LIOParams(..),
+        LIO(..),
+        invokeLIO,
+        runLIO,
+        lioParams,
+        rootLIOParams,
+        parentLIOParams,
+        upSideLIOParams,
+        downSideLIOParams,
+        shiftLIOParams,
+        latticeTimeIndex,
+        latticeMemberIndex,
+        latticeTime,
+        latticeTimeStep,
+        latticePoint) where
+
+import Data.IORef
+import Data.Maybe
+
+import Control.Applicative
+import Control.Monad
+import Control.Monad.Trans
+import Control.Monad.Fix
+import Control.Exception (throw, catch, finally)
+
+import Simulation.Aivika.Trans.Exception
+import Simulation.Aivika.Trans.Internal.Types
+import Simulation.Aivika.Trans.Parameter
+
+-- | The 'LIO' computation that can be run as nested one on the lattice node.
+newtype LIO a = LIO { unLIO :: LIOParams -> IO a
+                      -- ^ Unwrap the computation.
+                    }
+
+-- | The parameters of the 'LIO' computation.
+data LIOParams =
+  LIOParams { lioTimeIndex :: !Int,
+              -- ^ The time index.
+              lioMemberIndex :: !Int
+              -- ^ The member index.
+            } deriving (Eq, Ord, Show)
+
+instance Monad LIO where
+
+  {-# INLINE return #-}
+  return = LIO . const . return
+
+  {-# INLINE (>>=) #-}
+  (LIO m) >>= k = LIO $ \ps ->
+    m ps >>= \a ->
+    let m' = unLIO (k a) in m' ps
+
+instance Applicative LIO where
+
+  {-# INLINE pure #-}
+  pure = return
+
+  {-# INLINE (<*>) #-}
+  (<*>) = ap
+
+instance Functor LIO where
+
+  {-# INLINE fmap #-}
+  fmap f (LIO m) = LIO $ fmap f . m 
+
+instance MonadIO LIO where
+
+  {-# INLINE liftIO #-}
+  liftIO = LIO . const . liftIO
+
+instance MonadFix LIO where
+
+  mfix f = 
+    LIO $ \ps ->
+    do { rec { a <- invokeLIO ps (f a) }; return a }
+
+instance MonadException LIO where
+
+  catchComp (LIO m) h = LIO $ \ps ->
+    catch (m ps) (\e -> unLIO (h e) ps)
+
+  finallyComp (LIO m1) (LIO m2) = LIO $ \ps ->
+    finally (m1 ps) (m2 ps)
+  
+  throwComp e = LIO $ \ps ->
+    throw e
+
+-- | Invoke the computation.
+invokeLIO :: LIOParams -> LIO a -> IO a
+{-# INLINE invokeLIO #-}
+invokeLIO ps (LIO m) = m ps
+
+-- | Run the 'LIO' computation using the integration times points as lattice nodes.
+runLIO :: LIO a -> IO a
+runLIO m = unLIO m rootLIOParams
+
+-- | Return the parameters of the computation.
+lioParams :: LIO LIOParams
+lioParams = LIO return
+
+-- | Return the root node parameters.
+rootLIOParams :: LIOParams
+rootLIOParams = LIOParams { lioTimeIndex = 0,
+                            lioMemberIndex = 0 }
+
+-- | Return the parent parameters.
+parentLIOParams :: LIOParams -> Maybe LIOParams
+parentLIOParams ps
+  | i == 0    = Nothing
+  | otherwise = Just $ ps { lioTimeIndex = i - 1, lioMemberIndex = max 0 (k - 1) }
+  where i = lioTimeIndex ps
+        k = lioMemberIndex ps
+
+-- | Return the next up side parameters.
+upSideLIOParams :: LIOParams -> LIOParams
+upSideLIOParams ps = ps { lioTimeIndex = 1 + i }
+  where i = lioTimeIndex ps
+
+-- | Return the next down side parameters.
+downSideLIOParams :: LIOParams -> LIOParams
+downSideLIOParams ps = ps { lioTimeIndex = 1 + i, lioMemberIndex = 1 + k }
+  where i = lioTimeIndex ps
+        k = lioMemberIndex ps
+
+-- | Return the derived LIO params with the specified shift in 'latticeTimeIndex' and
+-- 'latticeMemberIndex' respectively, where the first parameter can be positive only.
+shiftLIOParams :: Int
+                  -- ^ a positive shift the lattice time index
+                  -> Int
+                  -- ^ a shift of the lattice member index
+                  -> LIOParams
+                  -- ^ the source parameters
+                  -> LIOParams
+shiftLIOParams di dk ps
+  | di <= 0   = error "The time index shift must be positive: shiftLIOParams"
+  | k' < 0    = error "The member index cannot be negative: shiftLIOParams"
+  | k' > i'   = error "The member index cannot be greater than the time index: shiftLIOParams"
+  | otherwise = ps { lioTimeIndex = i', lioMemberIndex = k' }
+  where i  = lioTimeIndex ps
+        i' = i + di
+        k  = lioMemberIndex ps
+        k' = k + dk
+
+-- | Return the lattice time index starting from 0. It corresponds to the integration time point.
+latticeTimeIndex :: LIO Int
+latticeTimeIndex = LIO $ return . lioTimeIndex
+
+-- | Return the lattice member index starting from 0. It is always less than or equaled to 'latticeTimeIndex'.
+latticeMemberIndex :: LIO Int
+latticeMemberIndex = LIO $ return . lioMemberIndex
+
+-- | Return the time for the current lattice node.
+latticeTime :: Parameter LIO Double
+latticeTime =
+  Parameter $ \r ->
+  LIO $ \ps ->
+  let sc = runSpecs r
+      i  = lioTimeIndex ps
+      t  = spcStartTime sc + (fromInteger $ toInteger i) * (spcDT sc)
+  in return t
+
+-- | Return the point in the corresponding lattice node.
+latticePoint :: Parameter LIO (Point LIO)
+latticePoint =
+  Parameter $ \r ->
+  do t <- invokeParameter r latticeTime
+     let sc = runSpecs r
+         t0 = spcStartTime sc
+         dt = spcDT sc
+         n  = fromIntegral $ floor ((t - t0) / dt)
+     return Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = t,
+                    pointIteration = n,
+                    pointPhase = -1 }
+
+-- | The time step used when constructing the lattice. Currently, it is equivalent to 'dt'.
+latticeTimeStep :: Parameter LIO Double
+latticeTimeStep = dt
diff --git a/Simulation/Aivika/Lattice/Internal/Ref.hs b/Simulation/Aivika/Lattice/Internal/Ref.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Internal/Ref.hs
@@ -0,0 +1,168 @@
+
+{-# LANGUAGE BangPatterns #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Internal.Ref
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- The implementation of mutable references.
+--
+module Simulation.Aivika.Lattice.Internal.Ref
+       (Ref,
+        newEmptyRef,
+        newEmptyRef0,
+        newRef,
+        newRef0,
+        readRef,
+        readRef0,
+        writeRef,
+        writeRef0,
+        modifyRef,
+        modifyRef0,
+        topRefDefined0,
+        defineTopRef0,
+        defineTopRef0_) where
+
+-- import Debug.Trace
+
+import Data.IORef
+import qualified Data.IntMap as M
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans.Internal.Types
+import Simulation.Aivika.Lattice.Internal.LIO
+
+-- | A reference map.
+type RefMap a = IORef (M.IntMap (IORef a))
+
+-- | A mutable reference.
+newtype Ref a = Ref { refMap :: RefMap a
+                      -- ^ the map of actual references
+                    }
+
+instance Eq (Ref a) where
+  r1 == r2 = (refMap r1) == (refMap r2)
+
+-- | Return the map index.
+lioMapIndex :: LIOParams -> Int
+lioMapIndex ps = ((i * (i + 1)) `div` 2) + k
+  where i = lioTimeIndex ps
+        k = lioMemberIndex ps
+
+-- | Create an empty reference.
+newEmptyRef :: Simulation LIO (Ref a)
+newEmptyRef = Simulation $ const newEmptyRef0
+
+-- | Create an empty reference.
+newEmptyRef0 :: LIO (Ref a)
+newEmptyRef0 =
+  LIO $ \ps ->
+  do rm <- newIORef M.empty
+     return Ref { refMap = rm }
+
+-- | Create a new reference.
+newRef :: a -> Simulation LIO (Ref a)
+newRef = Simulation . const . newRef0
+
+-- | Create a new reference.
+newRef0 :: a -> LIO (Ref a)
+newRef0 a =
+  LIO $ \ps ->
+  do r  <- invokeLIO ps newEmptyRef0
+     ra <- newIORef a
+     let !i = lioMapIndex ps
+     writeIORef (refMap r) $
+       M.insert i ra M.empty
+     return r
+     
+-- | Read the value of a reference.
+readRef :: Ref a -> Event LIO a
+readRef = Event . const . readRef0
+     
+-- | Read the value of a reference.
+readRef0 :: Ref a -> LIO a
+readRef0 r =
+  LIO $ \ps ->
+  do m <- readIORef (refMap r)
+     let loop ps =
+           case M.lookup (lioMapIndex ps) m of
+             Just ra -> readIORef ra
+             Nothing ->
+               case parentLIOParams ps of
+                 Just ps' -> loop ps'
+                 Nothing  -> error "Cannot find lattice node: readRef0"
+     loop ps
+
+-- | Write a new value into the reference.
+writeRef :: Ref a -> a -> Event LIO ()
+writeRef r a = Event $ const $ writeRef0 r a 
+
+-- | Write a new value into the reference.
+writeRef0 :: Ref a -> a -> LIO ()
+writeRef0 r a =
+  LIO $ \ps ->
+  do m <- readIORef (refMap r)
+     let !i = lioMapIndex ps
+     case M.lookup i m of
+       Just ra -> a `seq` writeIORef ra a
+       Nothing ->
+         do ra <- a `seq` newIORef a
+            modifyIORef (refMap r) $ M.insert i ra
+
+-- | Mutate the contents of the reference.
+modifyRef :: Ref a -> (a -> a) -> Event LIO ()
+modifyRef r f = Event $ const $ modifyRef0 r f
+
+-- | Mutate the contents of the reference.
+modifyRef0 :: Ref a -> (a -> a) -> LIO ()
+modifyRef0 r f =
+  LIO $ \ps ->
+  do m <- readIORef (refMap r)
+     let !i = lioMapIndex ps
+     case M.lookup i m of
+       Just ra ->
+         do a <- readIORef ra
+            let b = f a
+            b `seq` writeIORef ra b
+       Nothing ->
+         do a <- invokeLIO ps $ readRef0 r
+            invokeLIO ps $ writeRef0 r (f a)
+
+-- | Whether the top reference value is defined.
+topRefDefined0 :: Ref a -> LIO Bool
+topRefDefined0 r =
+  LIO $ \ps ->
+  do m <- readIORef (refMap r)
+     let !i = lioMapIndex ps
+     case M.lookup i m of
+       Just ra -> return True
+       Nothing -> return False
+
+-- | Define the top reference value.
+defineTopRef0 :: Ref a -> LIO a
+defineTopRef0 r =
+  LIO $ \ps ->
+  do m <- readIORef (refMap r)
+     let !i = lioMapIndex ps
+     case M.lookup i m of
+       Just ra -> readIORef ra
+       Nothing ->
+         case parentLIOParams ps of
+           Nothing  -> error "Cannot find parent: defineTopRef0"
+           Just ps' ->
+             do a  <- invokeLIO ps' $ defineTopRef0 r
+                ra <- a `seq` newIORef a
+                modifyIORef (refMap r) $ M.insert i ra
+                return a
+
+-- | Ensure that the top reference value is defined.
+defineTopRef0_ :: Ref a -> LIO ()
+defineTopRef0_ r =
+  do a <- defineTopRef0 r
+     return ()
diff --git a/Simulation/Aivika/Lattice/LIO.hs b/Simulation/Aivika/Lattice/LIO.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/LIO.hs
@@ -0,0 +1,40 @@
+
+-- |
+-- Module     : Simulation.Aivika.Branch.LIO
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- This module defines 'LIO' as an instance of the 'MonadDES' and 'EventIOQueueing' type classes.
+--
+module Simulation.Aivika.Lattice.LIO
+       (LIO,
+        runLIO,
+        latticeTimeIndex,
+        latticeMemberIndex) where
+
+import Simulation.Aivika.Trans.Comp
+import Simulation.Aivika.Trans.DES
+import Simulation.Aivika.Trans.Exception
+import Simulation.Aivika.Trans.Generator
+import Simulation.Aivika.Trans.Event
+import Simulation.Aivika.Trans.Ref.Base
+import Simulation.Aivika.Trans.QueueStrategy
+
+import Simulation.Aivika.Lattice.Internal.LIO
+import Simulation.Aivika.Lattice.Event
+import Simulation.Aivika.Lattice.Generator
+import Simulation.Aivika.Lattice.Ref.Base
+import Simulation.Aivika.Lattice.QueueStrategy
+
+instance MonadDES LIO
+
+instance MonadComp LIO
+
+-- | An implementation of the 'EventIOQueueing' type class.
+instance EventIOQueueing LIO where
+
+  enqueueEventIO = enqueueEvent
+  
diff --git a/Simulation/Aivika/Lattice/QueueStrategy.hs b/Simulation/Aivika/Lattice/QueueStrategy.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/QueueStrategy.hs
@@ -0,0 +1,68 @@
+
+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.QueueStrategy
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- This module defines queue strategies 'FCFS' and 'LCFS' for the 'LIO' computation.
+--
+module Simulation.Aivika.Lattice.QueueStrategy () where
+
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import qualified Simulation.Aivika.Trans.DoubleLinkedList as LL
+
+import Simulation.Aivika.Lattice.Internal.LIO
+import Simulation.Aivika.Lattice.Ref.Base
+
+-- | An implementation of the 'FCFS' queue strategy.
+instance QueueStrategy LIO FCFS where
+
+  -- | A queue used by the 'FCFS' strategy.
+  newtype StrategyQueue LIO FCFS a = FCFSQueueLIO (LL.DoubleLinkedList LIO a)
+
+  newStrategyQueue s = fmap FCFSQueueLIO LL.newList
+
+  strategyQueueNull (FCFSQueueLIO q) = LL.listNull q
+
+-- | An implementation of the 'FCFS' queue strategy.
+instance DequeueStrategy LIO FCFS where
+
+  strategyDequeue (FCFSQueueLIO q) =
+    do i <- LL.listFirst q
+       LL.listRemoveFirst q
+       return i
+
+-- | An implementation of the 'FCFS' queue strategy.
+instance EnqueueStrategy LIO FCFS where
+
+  strategyEnqueue (FCFSQueueLIO q) i = LL.listAddLast q i
+
+-- | An implementation of the 'LCFS' queue strategy.
+instance QueueStrategy LIO LCFS where
+
+  -- | A queue used by the 'LCFS' strategy.
+  newtype StrategyQueue LIO LCFS a = LCFSQueueLIO (LL.DoubleLinkedList LIO a)
+
+  newStrategyQueue s = fmap LCFSQueueLIO LL.newList
+       
+  strategyQueueNull (LCFSQueueLIO q) = LL.listNull q
+
+-- | An implementation of the 'LCFS' queue strategy.
+instance DequeueStrategy LIO LCFS where
+
+  strategyDequeue (LCFSQueueLIO q) =
+    do i <- LL.listFirst q
+       LL.listRemoveFirst q
+       return i
+
+-- | An implementation of the 'LCFS' queue strategy.
+instance EnqueueStrategy LIO LCFS where
+
+  strategyEnqueue (LCFSQueueLIO q) i = LL.listInsertFirst q i
diff --git a/Simulation/Aivika/Lattice/Ref/Base.hs b/Simulation/Aivika/Lattice/Ref/Base.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Lattice/Ref/Base.hs
@@ -0,0 +1,59 @@
+
+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}
+
+-- |
+-- Module     : Simulation.Aivika.Lattice.Ref.Base
+-- Copyright  : Copyright (c) 2016, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.10.3
+--
+-- Here is an implementation of mutable references, where
+-- 'LIO' is an instance of 'MonadRef' and 'MonadRef0'.
+--
+module Simulation.Aivika.Lattice.Ref.Base () where
+
+import Simulation.Aivika.Trans.Internal.Types
+import Simulation.Aivika.Trans.Comp
+import Simulation.Aivika.Trans.Simulation
+import Simulation.Aivika.Trans.Ref.Base
+import Simulation.Aivika.Trans.Observable
+
+import Simulation.Aivika.Lattice.Internal.LIO
+import qualified Simulation.Aivika.Lattice.Internal.Ref as R
+import Simulation.Aivika.Lattice.Internal.Estimate
+import Simulation.Aivika.Lattice.Internal.Event
+
+-- | The implementation of mutable references.
+instance MonadRef LIO where
+
+  -- | The mutable reference.
+  newtype Ref LIO a = Ref { refValue :: R.Ref a }
+
+  {-# INLINE newRef #-}
+  newRef = fmap Ref . R.newRef 
+
+  {-# INLINE readRef #-}
+  readRef (Ref r) = R.readRef r
+
+  {-# INLINE writeRef #-}
+  writeRef (Ref r) = R.writeRef r
+
+  {-# INLINE modifyRef #-}
+  modifyRef (Ref r) = R.modifyRef r
+
+  {-# INLINE equalRef #-}
+  equalRef (Ref r1) (Ref r2) = (r1 == r2)
+
+-- | A subtype of mutable references that can be created under more weak conditions.
+instance MonadRef0 LIO where
+
+  {-# INLINE newRef0 #-}
+  newRef0 = fmap Ref . R.newRef0
+
+-- | An instance of the specified type class.
+instance Observable (Ref LIO) (Estimate LIO) where
+
+  {-# INLINE readObservable #-}
+  readObservable (Ref r) = estimateRef r
diff --git a/aivika-lattice.cabal b/aivika-lattice.cabal
new file mode 100644
--- /dev/null
+++ b/aivika-lattice.cabal
@@ -0,0 +1,62 @@
+name:            aivika-lattice
+version:         0.1
+synopsis:        Nested discrete event simulation module for the Aivika library using lattice
+description:
+    This experimental package extends the aivika-transformers [1] library with facilities for 
+    running nested discrete event simulations within lattice nodes.
+    .
+    \[1] <http://hackage.haskell.org/package/aivika-transformers>
+    .
+category:        Simulation
+license:         BSD3
+license-file:    LICENSE
+copyright:       (c) 2016. David Sorokin <david.sorokin@gmail.com>
+author:          David Sorokin
+maintainer:      David Sorokin <david.sorokin@gmail.com>
+homepage:        http://www.aivikasoft.com/en/products/aivika.html
+cabal-version:   >= 1.6
+build-type:      Simple
+tested-with:     GHC == 7.10.3
+
+extra-source-files:  CHANGELOG.md
+                     tests/MachRep1.hs
+                     tests/TraversingLattice1.hs
+                     tests/TraversingLattice2.hs
+                     tests/TraversingLattice3.hs
+                     tests/TraversingLattice4.hs
+                     tests/TraversingMachRep1.hs
+                     tests/EstimatingMachRep1.hs
+
+library
+
+    exposed-modules: Simulation.Aivika.Lattice
+                     Simulation.Aivika.Lattice.Estimate
+                     Simulation.Aivika.Lattice.Event
+                     Simulation.Aivika.Lattice.Generator
+                     Simulation.Aivika.Lattice.LIO
+                     Simulation.Aivika.Lattice.QueueStrategy
+                     Simulation.Aivika.Lattice.Ref.Base
+
+    other-modules:   Simulation.Aivika.Lattice.Internal.Event
+                     Simulation.Aivika.Lattice.Internal.Estimate
+                     Simulation.Aivika.Lattice.Internal.LIO
+                     Simulation.Aivika.Lattice.Internal.Ref
+
+    build-depends:   base >= 3 && < 6,
+                     mtl >= 2.1.1,
+                     containers >= 0.4.0.0,
+                     random >= 1.0.0.3,
+                     aivika >= 4.5,
+                     aivika-transformers >= 4.5
+
+    extensions:      TypeFamilies,
+                     MultiParamTypeClasses,
+                     FlexibleInstances,
+                     BangPatterns
+
+    ghc-options:     -O2
+
+source-repository head
+
+    type:     git
+    location: https://github.com/dsorokin/aivika-lattice
diff --git a/tests/EstimatingMachRep1.hs b/tests/EstimatingMachRep1.hs
new file mode 100644
--- /dev/null
+++ b/tests/EstimatingMachRep1.hs
@@ -0,0 +1,76 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 100.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO Double
+model =
+  do totalUpTime <- newRef 0.0
+     
+     let machine =
+           do upTime <-
+                liftParameter $
+                randomExponential meanUpTime
+              --
+              -- r <- liftSimulation $ newRef 10
+              --
+              holdProcess upTime
+              liftEvent $ 
+                modifyRef totalUpTime (+ upTime)
+              repairTime <-
+                liftParameter $
+                randomExponential meanRepairTime
+              holdProcess repairTime
+              machine
+
+     runProcessInStartTime machine
+     runProcessInStartTime machine
+
+     t2 <- liftParameter stoptime
+
+     let leaf =
+           do x <- readObservable totalUpTime
+              t <- estimateTime
+              let a = x / (2 * t)
+              return a
+
+     let reduce a1 a2 =
+           do let a = (a1 + a2) / 2
+              return a
+
+     r <- foldEstimate reduce leaf
+
+     runEstimateInStartTime r
+
+main :: IO ()
+main =
+  do a <- runLIO $
+          runSimulation model specs
+     print a
diff --git a/tests/MachRep1.hs b/tests/MachRep1.hs
new file mode 100644
--- /dev/null
+++ b/tests/MachRep1.hs
@@ -0,0 +1,73 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 100.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO (Results LIO)
+model =
+  do totalUpTime <- newRef 0.0
+     
+     let machine =
+           do upTime <-
+                liftParameter $
+                randomExponential meanUpTime
+              --
+              -- r <- liftSimulation $ newRef 10
+              --
+              holdProcess upTime
+              liftEvent $ 
+                modifyRef totalUpTime (+ upTime)
+              repairTime <-
+                liftParameter $
+                randomExponential meanRepairTime
+              holdProcess repairTime
+              machine
+
+     runProcessInStartTime machine
+     runProcessInStartTime machine
+
+     let upTimeProp =
+           do x <- readRef totalUpTime
+              t <- liftDynamics time
+              return $ x / (2 * t)
+
+     return $
+       results
+       [resultSource
+        "upTimeProp"
+        "The long-run proportion of up time (~ 0.66)"
+        upTimeProp]
+
+main :: IO ()
+main =
+  runLIO $
+  printSimulationResultsInStopTime
+  printResultSourceInEnglish
+  model specs
diff --git a/tests/TraversingLattice1.hs b/tests/TraversingLattice1.hs
new file mode 100644
--- /dev/null
+++ b/tests/TraversingLattice1.hs
@@ -0,0 +1,40 @@
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 100.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO ()
+model =
+  do let showLatticeNode :: Event LIO ()
+         showLatticeNode =
+           do t <- liftDynamics time
+              i <- liftComp latticeTimeIndex
+              k <- liftComp latticeMemberIndex
+              liftIO $
+                do putStr $ "t = " ++ show t
+                   putStr $ ", time index = " ++ show i
+                   putStr $ ", member index = " ++ show k
+                   putStrLn ""
+                   
+     runEventInStartTime $
+       enqueueEventWithIntegTimes
+       showLatticeNode
+
+     runEventInStopTime $
+       return ()
+
+main :: IO ()
+main =
+  runLIO $
+  runSimulation model specs
diff --git a/tests/TraversingLattice2.hs b/tests/TraversingLattice2.hs
new file mode 100644
--- /dev/null
+++ b/tests/TraversingLattice2.hs
@@ -0,0 +1,52 @@
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 400.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO ()
+model =
+  do let showLatticeNode :: String -> Event LIO ()
+         showLatticeNode action =
+           do t <- liftDynamics time
+              i <- liftComp latticeTimeIndex
+              k <- liftComp latticeMemberIndex
+              liftIO $
+                do putStr action
+                   putStr $ ": t = " ++ show t
+                   putStr $ ", time index = " ++ show i
+                   putStr $ ", member index = " ++ show k
+                   putStrLn ""
+                   
+     runEventInStartTime $
+       enqueueEventWithIntegTimes $
+       showLatticeNode "enqueue"
+
+     let reduce a b =
+           traceEstimate "reduce" $
+           return ()
+
+     let leaf =
+           traceEstimate "leaf" $
+           return ()
+
+     foldEstimate reduce leaf
+       >>= runEstimateInStartTime
+
+     runEventInStopTime $
+       return ()
+
+main :: IO ()
+main =
+  runLIO $
+  runSimulation model specs
diff --git a/tests/TraversingLattice3.hs b/tests/TraversingLattice3.hs
new file mode 100644
--- /dev/null
+++ b/tests/TraversingLattice3.hs
@@ -0,0 +1,52 @@
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 400.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO ()
+model =
+  do let showLatticeNode :: String -> Event LIO ()
+         showLatticeNode action =
+           do t <- liftDynamics time
+              i <- liftComp latticeTimeIndex
+              k <- liftComp latticeMemberIndex
+              liftIO $
+                do putStr action
+                   putStr $ ": t = " ++ show t
+                   putStr $ ", time index = " ++ show i
+                   putStr $ ", member index = " ++ show k
+                   putStrLn ""
+                   
+     runEventInStartTime $
+       enqueueEventWithIntegTimes $
+       showLatticeNode "enqueue"
+
+     let reduce a b =
+           traceEstimate "reduce" $
+           return ()
+
+     let leaf =
+           traceEstimate "leaf" $
+           return ()
+
+     foldEstimate reduce leaf
+       >>= runEstimateInStartTime
+
+     runEventInStopTime $
+       showLatticeNode "stop"
+
+main :: IO ()
+main =
+  runLIO $
+  runSimulation model specs
diff --git a/tests/TraversingLattice4.hs b/tests/TraversingLattice4.hs
new file mode 100644
--- /dev/null
+++ b/tests/TraversingLattice4.hs
@@ -0,0 +1,63 @@
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 400.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO ()
+model =
+  do let showLatticeNode :: String -> Event LIO ()
+         showLatticeNode action =
+           do t <- liftDynamics time
+              i <- liftComp latticeTimeIndex
+              k <- liftComp latticeMemberIndex
+              liftIO $
+                do putStr action
+                   putStr $ ": t = " ++ show t
+                   putStr $ ", time index = " ++ show i
+                   putStr $ ", member index = " ++ show k
+                   putStrLn ""
+
+     r <- newRef 0
+                   
+     runEventInStartTime $
+       enqueueEventWithIntegTimes $
+       do x <- liftParameter $ randomUniform 0 1
+          writeRef r x
+          showLatticeNode ("enqueue (x = " ++ show x ++ ")") 
+
+     let reduce :: Double -> Double -> Estimate LIO Double
+         reduce a b =
+           do let x = (a + b) / 2
+              traceEstimate ("reduce (x = " ++ show x ++ ")") $
+                return x
+
+     let leaf =
+           do x <- readObservable r
+              traceEstimate ("leaf (x = " ++ show x ++ ")") $
+                return x
+
+     m <- foldEstimate reduce leaf
+
+     runEstimateInStartTime $
+       do x <- m
+          traceEstimate ("result (x = " ++ show x ++ ")") $
+            return ()
+
+     runEventInStopTime $
+       showLatticeNode "stop"
+
+main :: IO ()
+main =
+  runLIO $
+  runSimulation model specs
diff --git a/tests/TraversingMachRep1.hs b/tests/TraversingMachRep1.hs
new file mode 100644
--- /dev/null
+++ b/tests/TraversingMachRep1.hs
@@ -0,0 +1,79 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Trans
+import Simulation.Aivika.Lattice
+
+meanUpTime = 1.0
+meanRepairTime = 0.5
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 100.0,
+                spcMethod = RungeKutta4,
+                spcGeneratorType = SimpleGenerator }
+        
+model :: Simulation LIO Double
+model =
+  do totalUpTime <- newRef 0.0
+     
+     let machine =
+           do upTime <-
+                liftParameter $
+                randomExponential meanUpTime
+              --
+              -- r <- liftSimulation $ newRef 10
+              --
+              holdProcess upTime
+              liftEvent $ 
+                modifyRef totalUpTime (+ upTime)
+              repairTime <-
+                liftParameter $
+                randomExponential meanRepairTime
+              holdProcess repairTime
+              machine
+
+     runProcessInStartTime machine
+     runProcessInStartTime machine
+
+     t2 <- liftParameter stoptime
+
+     let leaf =
+           do x <- readObservable totalUpTime
+              t <- estimateTime
+              let a = x / (2 * t)
+              traceEstimate ("leaf (" ++ show a ++ ")") $
+                return a
+
+     let reduce a1 a2 =
+           do let a = (a1 + a2) / 2
+              traceEstimate ("result (" ++ show a ++ ")") $ 
+                return a
+
+     r <- foldEstimate reduce leaf
+
+     runEstimateInStartTime $
+       traceEstimate "launch" r
+
+main :: IO ()
+main =
+  do a <- runLIO $
+          runSimulation model specs
+     print a