diff --git a/src/Data/Dim.hs b/src/Data/Dim.hs
--- a/src/Data/Dim.hs
+++ b/src/Data/Dim.hs
@@ -1,4 +1,5 @@
 {-# language GADTs, TypeOperators #-}
+{-# language TypeFamilies #-}
 module Data.Dim where
 
 import Data.Vector.Unboxed as V
@@ -6,10 +7,13 @@
 
 -- * Dimension metadata
 
+data family Dim a
+data instance Dim (Dd i)
+data instance Dim (Sd i)
+
 -- | To define a /dense/ dimension we only need the dimensionality parameter
 newtype Dd i = Dd { dDim :: i } deriving (Eq, Show)
 
-
 -- | To define a /sparse/ dimension we need a cumulative array, an index array and a dimensionality parameter
 data Sd i = Sd {
       -- | Cumulative array (# nonzero entries per degree of freedom). Not all storage formats (e.g. COO for rank-2 tensors) need this information.
@@ -24,7 +28,8 @@
 --
 -- Example: the CSR format is /dense/ in the first index (rows) and /sparse/ in the second index (columns)
 
--- newtype DMD i = DMD (Either (DMDDense i) (DMDSparse i)) deriving (Eq, Show)
+dim :: Either (Dd c) (Sd c) -> c
+dim = either dDim sDim
 
 
   
diff --git a/src/Data/Shape.hs b/src/Data/Shape.hs
--- a/src/Data/Shape.hs
+++ b/src/Data/Shape.hs
@@ -1,4 +1,5 @@
 {-# language GADTs, TypeOperators #-}
+{-# language TypeFamilies #-}
 module Data.Shape where
 
 import Data.Monoid
@@ -12,12 +13,24 @@
 import qualified Data.Dim as Dim
    
 
+-- | A class for data that have a shape.
+class Shape t where
+  type ShapeT t :: *
+  shape :: t -> ShapeT t
+  shRank :: t -> Int
+  shDim :: t -> [Int]
+
+
+
 data Z
 data sh :# e -- dense
 data sh :. e -- sparse
 
 -- | A statically-typed tensor shape parameter that supports both sparse and dense dimensions.
 -- Dimensions are indexed with 'Int32' indices, which should be enough for most applications.
+--
+-- Note: in this formulation, only the /rank/ (i.e. the number of dimensions) is known at compile time. The actual dimensionality is only known at runtime.
+-- To address this, see Data.Shape.Static
 data Sh sh where 
   Z :: Sh Z
   -- | Constructor for a dense dimension  
diff --git a/src/Data/Shape/Dynamic.hs b/src/Data/Shape/Dynamic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Shape/Dynamic.hs
@@ -0,0 +1,23 @@
+module Data.Shape.Dynamic where
+
+import Data.Foldable (foldl')
+
+import qualified Data.Dim as Dim
+
+
+newtype ShD i =
+    ShD {unShD :: [Either (Dim.Dd i) (Dim.Sd i)] }
+      deriving (Eq, Show)
+
+mkDenseShD :: Foldable t => t i -> ShD i
+mkDenseShD xss = ShD $ foldl' (\d s -> Left (Dim.Dd s) : d) [] xss
+
+rank :: ShD i -> Int
+rank = length . unShD 
+
+-- dim :: ShD i -> [i]
+dim :: (Num b, Integral a) => ShD a -> [b]
+dim sh = fromIntegral <$> foldl' (\d s -> Dim.dim s : d) [] (unShD sh)
+
+
+
diff --git a/src/Data/Shape/Static.hs b/src/Data/Shape/Static.hs
--- a/src/Data/Shape/Static.hs
+++ b/src/Data/Shape/Static.hs
@@ -1,11 +1,138 @@
 {-# language GADTs, TypeOperators, DataKinds, KindSignatures #-}
+{-# language FlexibleInstances, RankNTypes #-}
+{-# language TypeFamilies #-}
+{-# language PolyKinds #-}
+{-# language MultiParamTypeClasses #-}
+-- {-# language TypeInType #-}
 module Data.Shape.Static where
 
-import           GHC.TypeLits (Nat)
+import Data.Monoid ((<>))
+import Data.Int (Int32)
 
--- from https://hackage.haskell.org/package/dimensions-0.3.2.0/docs/src/Numeric-Dimensions-Idx.html#Idx
-data Idx (ds :: [Nat]) where
-   -- | Zero-rank dimensionality - scalar
-   Zi :: Idx '[]
-   -- | List-like concatenation of indices
-   (:!) :: {-# UNPACK #-} !Int -> !(Idx ds) -> Idx (d ': ds)
+import GHC.TypeLits -- (Nat, natVal, KnownNat(..))
+import Data.Proxy (Proxy(..))
+import Unsafe.Coerce (unsafeCoerce)
+
+import GHC.Generics (Generic)
+
+import qualified Data.Vector.Unboxed as VU
+
+import qualified Data.Dim as Dim
+
+
+
+
+ 
+data sh :# e -- dense
+data sh :. e -- sparse
+infixr 5 :#
+infixr 5 :.
+
+
+data D (n :: Nat)
+data S (n :: Nat)
+data family Dimz d
+data instance Dimz (D n)
+data instance Dimz (S n)
+
+-- | A shape type with statically typed dimensions  
+data Sh sh where
+  Z :: Sh '[]
+  -- | Constructor for a dense dimension
+  D :: KnownNat n => Dim.Dd Int32 -> Sh sh -> Sh (D n ': sh)
+  -- | Constructor for a sparse dimension
+  S :: KnownNat n => Dim.Sd Int32 -> Sh sh -> Sh (S n ': sh)
+
+-- data Sized (n :: [Dimz Nat]) t = Sized t   -- this doesn't work without TypeInType
+
+-- | 
+data Sized n c t = Sized (Sh n) (c t)
+
+t0 :: Sized '[D 3] [] Int
+t0 = Sized (Dim.Dd 3 `D` Z) [1,2,3]
+    
+  
+  
+instance Show (Sh sh) where
+  show Z = ""
+  show (D (Dim.Dd m) sh) = unwords [show m, show sh]
+  show (S (Dim.Sd _ ix n) sh) = showSparse ix n <> show sh where
+    showSparse ixx nn = show (VU.length ixx, nn)
+
+shToList :: Sh ds -> [Int32]
+shToList Z = []
+shToList (x `D` xs) = Dim.dDim x : shToList xs
+shToList (x `S` xs) = Dim.sDim x : shToList xs
+
+
+
+
+-- * `dimensions`
+
+-- data Sh (ds :: [Nat]) where
+--   Z :: Sh '[]
+--   D :: {-# UNPACK#-} !Int -> !(Sh ds) -> Sh (d ': ds)
+
+-- -- from https://hackage.haskell.org/package/dimensions-0.3.2.0/docs/src/Numeric-Dimensions-Idx.html#Idx
+-- data Idx (ds :: [Nat]) where
+--    -- | Zero-rank dimensionality - scalar
+--    Z :: Idx '[]
+--    -- | List-like concatenation of indices
+--    (:!) :: {-# UNPACK #-} !Int -> !(Idx ds) -> Idx (d ': ds)
+
+-- infixr 5 :!
+
+-- idxToList :: Idx ds -> [Int]
+-- idxToList Z = []
+-- idxToList (x :! xs) = x : idxToList xs
+
+-- -- | UNSAFE coerce 
+-- idxFromList :: [Int] -> Idx ds
+-- idxFromList xss = unsafeCoerce $ go xss
+--   where
+--     go [] = unsafeCoerce Z
+--     go (x:xs) = x :! unsafeCoerce (idxFromList xs)
+
+
+
+-- instance Show (Idx ds) where
+--     show Z  = "Idx Ø"
+--     show xs = "Idx" ++ foldr (\i s -> " " ++ show i ++ s) "" (idxToList xs)
+
+-- instance Eq (Idx ds) where
+--     Z == Z = True
+--     (a:!as) == (b:!bs) = a == b && as == bs
+--     Z /= Z = False
+--     (a:!as) /= (b:!bs) = a /= b || as /= bs
+
+
+-- -- | With this instance we can slightly reduce indexing expressions
+-- --   e.g. x ! (1 :! 2 :! 4) == x ! (1 :! 2 :! 4 :! Z)
+-- instance Num (Idx '[n]) where
+--     (a:!Z) + (b:!Z) = (a+b) :! Z
+--     (a:!Z) - (b:!Z) = (a-b) :! Z
+--     (a:!Z) * (b:!Z) = (a*b) :! Z
+--     signum (a:!Z)   = signum a :! Z
+--     abs (a:!Z)      = abs a :! Z
+--     fromInteger i   = fromInteger i :! Z
+
+-- instance Ord (Idx ds) where
+--     compare Z Z             = EQ
+--     compare (a:!as) (b:!bs) = compare as bs `mappend` compare a b
+
+
+
+
+-- instance Dimensions ds => Bounded (Idx ds) where
+--     maxBound = f (dim @ds)
+--       where
+--         f :: forall ns . Dim ns -> Idx ns
+--         f D                     = Z
+--         f ((Dn :: Dim n) :* ds) = dimVal' @n :! f ds
+--     {-# INLINE maxBound #-}
+--     minBound = f (dim @ds)
+--       where
+--         f :: forall (ns :: [Nat]) . Dim ns -> Idx ns
+--         f D          = Z
+--         f (Dn :* ds) = 1 :! f ds
+--     {-# INLINE minBound #-}    
diff --git a/src/Data/Tensor.hs b/src/Data/Tensor.hs
--- a/src/Data/Tensor.hs
+++ b/src/Data/Tensor.hs
@@ -2,7 +2,8 @@
 {-# language DeriveFunctor #-}
 {-# language TypeOperators #-}
 {-# language PackageImports #-}
-
+{-# language TypeFamilies #-}
+{-# language FlexibleInstances #-}
 module Data.Tensor (
   -- * Tensor type
   Tensor(..),
@@ -20,15 +21,32 @@
 import Control.Applicative
 
 
-import qualified Data.Shape as Shape (dim, rank)
+import qualified Data.Shape as Shape (Shape(..), dim, rank)
 import Data.Shape (Sh(..), 
                    Z,
-                   D1, D2, CSR, COO, mkD2, mkCSR, mkCOO) 
+                   D1, D2, CSR, COO, mkD2, mkCSR, mkCOO)
+import qualified Data.Shape.Dynamic as ShDyn       
 import qualified Data.Dim as Dim
 
 
 
--- | The 'Tensor' type. Tensor data entries are stored as one single array
+-- | A tensor type with dimensions only known at runtime
+data Tenzor i a = Tenzor {tzShape :: ShDyn.ShD i, tzData :: V.Vector a }
+
+instance Integral i => Shape.Shape (Tenzor i a) where
+  type ShapeT (Tenzor i a) = ShDyn.ShD i
+  shape = tzShape 
+  shRank = ShDyn.rank . Shape.shape
+  shDim = ShDyn.dim . Shape.shape
+
+instance Shape.Shape (Tensor (Sh i) a) where
+  type ShapeT (Tensor (Sh i) a) = Sh i
+  shape = tshape 
+  shRank = rank
+  shDim = dim
+
+
+-- | The 'Tensor' type with statically known shape. Tensor data entries are stored as one single array
 data Tensor i a where
   Tensor :: Sh i -> V.Vector a -> Tensor (Sh i) a 
 
diff --git a/src/Data/Tensor/Compiler.hs b/src/Data/Tensor/Compiler.hs
--- a/src/Data/Tensor/Compiler.hs
+++ b/src/Data/Tensor/Compiler.hs
@@ -2,76 +2,78 @@
 {-# language PackageImports #-}
 
 module Data.Tensor.Compiler (
-    contract
-    -- * Tensor types
-  , Tensor(..), Sh(..), Dd(..), Sd(..)
+  --   contract
+  --   -- * Tensor types
+  -- , Tensor(..), Sh(..), Dd(..), Sd(..)
     -- * Syntax
-  , Phoas, eval, var, let_, let2_
+    Phoas, eval, var, let_, let2_
     -- * Exceptions
   , CException (..)
   )where
 
 import Data.Typeable
-import "exceptions" Control.Monad.Catch (MonadThrow(..), throwM, MonadCatch(..), catch)
+-- import "exceptions" Control.Monad.Catch (MonadThrow(..), throwM, MonadCatch(..), catch)
 import Control.Exception (Exception(..))
 
-import Control.Applicative (liftA2, (<|>))
+-- import Control.Applicative (liftA2, (<|>))
 
-import Data.Tensor (Tensor(..), Sh(..), Dd(..), Sd(..), tshape, tdata, nnz, rank, dim)
-import Data.Tensor.Compiler.PHOAS (Phoas(..), let_, let2_, var, lift1, lift2, eval)
+import Data.Tensor -- (Tensor(..), Sh(..), Dd(..), Sd(..), tshape, tdata, nnz, rank, dim)
+import Data.Tensor.Compiler.PHOAS -- (Phoas(..), let_, let2_, var, lift1, lift2, eval)
 
 
-{- |
-IN: Tensor reduction syntax (Einstein notation)
 
-OUT: stride program (how to read/write memory)
+-- IN: Tensor reduction syntax (Einstein notation)
 
+-- OUT: stride program (how to read/write memory)
 
-taco compiles a tensor expression (e.g. C = A_{ijk}B_{k} ) into a series of nested loops.
 
-dimensions : can be either dense or sparse
+-- taco compiles a tensor expression (e.g. C = A_{ijk}B_{k} ) into a series of nested loops.
 
-internally, tensor data is stored in /dense/ vectors
+-- dimensions : can be either dense or sparse
 
-"contract A_{ijk}B_{k} over the third index"
+-- internally, tensor data is stored in /dense/ vectors
 
--}
+-- "contract A_{ijk}B_{k} over the third index"
 
 
 
 
-mkVar :: MonadThrow m => [Int] -> Tensor i a -> m (Phoas (Tensor i a))
-mkVar ixs0 t = do
-  ixs <- mkIxs ixs0 (rank t)
-  return $ var t
+-- | 
 
+-- mkVar :: MonadThrow m => [Int] -> Tensor i a -> m (Phoas (Tensor i a))
+-- mkVar ixs0 t = do
+--   ixs <- mkIxs ixs0 (rank t)
+--   return $ var t
 
-mkIxs :: MonadThrow m => [Int] -> Int -> m [Int]
-mkIxs ixs mm = go ixs []
-  where
-    go [] acc = pure acc
-    go (i:is) acc | i < 0 =
-                    throwM $ IncompatIx "Index must be non-negative"
-                  | i > mm - 1 =
-                    throwM $ IncompatIx $ unwords ["Index must be smaller than", show mm]
-                  | otherwise = go is (i : acc)
 
+-- mkIxs :: MonadThrow m => [Int] -> Int -> m [Int]
+-- mkIxs ixs mm = go ixs []
+--   where
+--     go [] acc = pure acc
+--     go (i:is) acc | i < 0 =
+--                     throwM $ IncompatIx "Index must be non-negative"
+--                   | i > mm - 1 =
+--                     throwM $ IncompatIx $ unwords ["Index must be smaller than", show mm]
+--                   | otherwise = go is (i : acc)
+                  
+
 -- | Tensor contraction
 --
 -- Inject two 'Tensor' constant into 'Var's, while ensuring that all the contraction indices are compatible with those of the tensors.
 --
 -- Throws a 'CException' if any index is nonnegative or too large for the shape of the given tensor.
-contract :: MonadThrow m =>
-                  [Int]           -- ^ Tensor contraction indices
-                  -> Tensor i a
-                  -> Tensor i b
-                  -> ([Int] -> Tensor i a -> Tensor i b -> Phoas c) -- ^ Contraction function
-                  -> m (Phoas c)
-contract ixs0 t1 t2 f = do
-  _ <- mkIxs ixs0 (rank t1)
-  ixs <- mkIxs ixs0 (rank t2)
-  pure $ let_ (var ixs) $ \ixs' ->
-    let2_ (var t1) (var t2) (f ixs')
+
+-- -- contract :: MonadThrow m =>
+-- --                   [Int]           -- ^ Tensor contraction indices
+-- --                   -> Tensor i1 a
+-- --                   -> Tensor i2 b
+-- --                   -> ([Int] -> Tensor i1 a -> Tensor i2 b -> Phoas c) -- ^ Contraction function
+-- --                   -> m (Phoas c)
+-- contract ixs0 t1 t2 f = do
+--   _ <- mkIxs ixs0 (rank t1)
+--   ixs <- mkIxs ixs0 (rank t2)
+--   pure $ let_ (var ixs) $ \ixs' ->
+--     let2_ (var t1) (var t2) (f ixs')
 
 
 -- | Exceptions
diff --git a/src/Data/Tensor/Compiler/MultiStage.hs b/src/Data/Tensor/Compiler/MultiStage.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Tensor/Compiler/MultiStage.hs
@@ -0,0 +1,66 @@
+{-# language GADTs #-}
+module Data.Tensor.Compiler.MultiStage where
+
+{-|
+* inner product of two vectors
+* matrix-vector action
+* matrix-matrix product
+-}
+
+-- | "User-facing" syntax
+data E1 a =
+    K1 a
+  | Dot a a -- (E1 a) (E1 a)
+  deriving (Eq, Show)
+
+-- | Evaluate user-facing syntax into internal one
+evalE1 :: E1 a -> E2 a
+evalE1 expr = case expr of
+  K1 x -> K2 x
+  Dot a b -> Fold Add (ZipWith Mul (K2 a) (K2 b))
+
+-- k1 :: a -> E1 a
+-- k1 = K1
+-- dot :: a -> a -> E1 a
+-- dot = Dot
+
+
+-- | Internal representation
+data E2 a =
+    K2 a
+  | ZipWith BinOp (E2 a) (E2 a)
+  | Fold BinOp (E2 a) deriving (Eq, Show)
+
+-- | Evaluate the internal representation into a concrete result. NB: the internal and returned type, [b], is constrained by the evaluation functions `zipWith`/`foldr`. In general this would be an opaque tensor data structure
+evalE2 :: Fractional b => b -> E2 [b] -> [b]
+evalE2 z expr = case expr of
+  K2 x -> x
+  ZipWith op v1 v2 -> zipWith (evalBinOp op) (evalE2 z v1) (evalE2 z v2)
+  Fold op v -> [foldr (evalBinOp op) z (evalE2 z v)]
+
+-- | E2 combinators
+k2 :: a -> E2 a
+k2 = K2
+addE2 :: a -> a -> E2 a
+addE2 a b = ZipWith Add (k2 a) (k2 b)
+mulE2 :: a -> a -> E2 a
+mulE2 a b = ZipWith Mul (k2 a) (k2 b)
+
+-- | Chained E1 -> E2 evaluation
+evalE :: Fractional b => b -> E1 [b] -> [b]
+evalE z = evalE2 z . evalE1
+
+
+
+-- data UnOp = Sqrt deriving (Eq, Show)
+data BinOp = Add | Sub | Mul | Div deriving (Eq, Show)
+evalBinOp :: Fractional a => BinOp -> a -> a -> a
+evalBinOp op = case op of
+  Add -> (+)
+  Sub -> (-)
+  Mul -> (*)
+  Div -> (/)
+
+-- v1, v2 :: E1 [Int]
+-- v1 = k [1..5]
+-- v2 = k [3..7]
diff --git a/src/Data/Tensor/Compiler/PHOAS.hs b/src/Data/Tensor/Compiler/PHOAS.hs
--- a/src/Data/Tensor/Compiler/PHOAS.hs
+++ b/src/Data/Tensor/Compiler/PHOAS.hs
@@ -84,42 +84,8 @@
 -- * A possible abstract syntax
 
 
-{-|
-* inner product of two vectors
-* matrix-vector action
-* matrix-matrix product
--}
 
-data Expr a =
-    Konst a
-  | Dot (Expr a) (Expr a)
-  deriving (Eq, Show)
 
-k :: a -> Expr a
-k = Konst
-
-
-
-
--- data UnOp = Sqrt deriving (Eq, Show)
-data BinOp = Add | Sub | Mul | Div deriving (Eq, Show)
-evalBinOp :: Fractional a => BinOp -> a -> a -> a
-evalBinOp op = case op of
-  Add -> (+)
-  Sub -> (-)
-  Mul -> (*)
-  Div -> (/)
-
-data EI a =
-    -- CW1 UnOp (Expr a)
-  CW2 BinOp (Expr a) (Expr a) deriving (Eq, Show)
-
-v1, v2 :: Expr [Int]
-v1 = k [1..5]
-v2 = k [3..7]
-
--- evalEI expr = case expr of
---   CW2 op x y ->
 
 
 
diff --git a/taco.cabal b/taco.cabal
--- a/taco.cabal
+++ b/taco.cabal
@@ -1,7 +1,7 @@
 name:                taco
-version:             0.2.0.0
+version:             0.3.0.0
 synopsis:            Tensor Algebra COmpiler
-description:         This library provides types and a compiler for tensor expressions.
+description:         Types and a compiler for tensor expressions.
 homepage:            https://github.com/ocramz/taco-hs#readme
 license:             BSD3
 license-file:        LICENSE
@@ -20,8 +20,10 @@
   other-modules:       Data.Tensor
                        Data.Shape
                        Data.Dim
+                       Data.Shape.Dynamic
                        Data.Shape.Static
                        Data.Tensor.Compiler.PHOAS
+                       Data.Tensor.Compiler.MultiStage
   build-depends:       base >= 4.7 && < 5
                      , exceptions
                      , mtl