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/bamstats.cabal b/bamstats.cabal
new file mode 100644
--- /dev/null
+++ b/bamstats.cabal
@@ -0,0 +1,21 @@
+Name:          bamstats
+Version:       0.1
+License:       GPL
+Cabal-Version: >= 1.6
+Build-Type:    Simple
+Category:      Bioinformatics
+Author:        Ketil Malde
+Maintainer:    Ketil Malde <ketil@malde.org>
+Synopsis:      A program to extract various information from BAM alignmnet files
+
+Library
+    Hs-Source-Dirs:  src
+    Exposed-Modules: Bio.SamTools.Classify
+    Ghc-Options:     -Wall
+
+Executable bam
+    Hs-Source-Dirs:  src
+    Main-Is:         Bam.hs
+    Build-Depends:   base >= 4 && < 5, samtools >= 0.2.1, cmdargs
+    Ghc-Options:     -Wall
+
diff --git a/src/Bam.hs b/src/Bam.hs
new file mode 100644
--- /dev/null
+++ b/src/Bam.hs
@@ -0,0 +1,40 @@
+-- | BAM - extract information from BAM files
+
+{-# Language DeriveDataTypeable #-}
+
+module Main where
+import Bio.SamTools.Classify
+import Bio.SamTools.Bam
+import System.Console.CmdArgs
+
+data Options = Stats { numrd :: Maybe Int, inputs :: [FilePath] }
+             | Hist { numrd :: Maybe Int, bins, maxdist :: Int, inputs :: [FilePath] }
+             | Dump { numrd :: Maybe Int, inputs :: [FilePath] }
+             deriving (Data,Typeable,Read,Show)
+
+clfy, hst, dump :: Options
+clfy = Stats { numrd = Nothing &= help "max number of reads (default: all)"
+                , inputs = []     &= args &= typ "BAM file(s)"
+                } &= help "Calculate statistics on insert sizes"
+hst  = Hist { numrd = Nothing &= help "max number of reads to include"
+            , bins  = 20      &= help "number of bins" &= typ "INT"
+            , maxdist = 1000  &= help "max insert size to count" &= typ "INT"
+            , inputs = []     &= args &= typ "BAM file(s)"
+            } &= help "Collect a histogram of insert sizes"
+dump = Dump { numrd = Just 100 &= help "max number of reads (default: 100)" 
+            , inputs = []     &= args &= typ "BAM file(s)"            
+            } &= help "Dump alignments in the different classes."
+       
+histgen o = H 0 [(b,0) | b <- enumFromThenTo size (2*size) (maxdist o)]
+  where size = (maxdist o `div` bins o)
+
+main :: IO ()
+main = do
+  o <- cmdArgs $ modes [clfy,hst,dump] 
+       &= help "Extract information from BAM files" 
+       &= program "bam" &= summary "bam v0.0, ©2011 Ketil Malde"
+  let geninp f = (case numrd o of Just x -> take x; Nothing -> id) `fmap` readBams f
+      genout = case o of Stats {} -> putStrLn . display . classify (cdef :: ClassStats)
+                         Hist {}  -> putStrLn . display . classify (histgen o)
+                         Dump {}  -> putStrLn . display . classify (cdef :: Collect)
+  mapM_ (\f -> genout =<< geninp f) $ inputs o
diff --git a/src/Bio/SamTools/Classify.hs b/src/Bio/SamTools/Classify.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/SamTools/Classify.hs
@@ -0,0 +1,132 @@
+
+-- | Stats reads in a BAM-file as innie outie rightie leftie,
+--   with stats for distances
+
+module Bio.SamTools.Classify where
+
+import Bio.SamTools.Bam
+import Data.List (foldl', intercalate)
+import Data.Maybe (isNothing, fromJust)
+import Text.Printf (printf)
+
+data Stats a = Class { total :: !Int, innies, outies, lefties, righties :: !a } deriving Show
+
+-- --------------------------------------------------
+-- The generic collection framework
+-- --------------------------------------------------
+
+-- | Pretty-print a 'Stats' value.
+display :: Insertable x => Stats x -> String
+display c = unlines $ map (intercalate "\t")
+  [  "Alignment" : dispheader (innies c)
+  ,  "innies   " : disp1 t (innies c)
+  ,  "outies   " : disp1 t (outies c)
+  ,  "lefties  " : disp1 t (lefties c)
+  ,  "righties " : disp1 t (righties c)
+  , ["Total reads: ", show t]]
+  where t = total c
+
+class Insertable x where
+  insert :: Bam1 -> x -> x
+  disp1  :: Int  -> x -> [String]
+  dispheader :: x -> [String]
+  cdef :: x  
+
+-- | Extract info from alignments
+classify :: Insertable x => x -> [Bam1] -> Stats x
+classify def = foldl' (class1 . bump) (Class 0 def def def def)
+
+-- | Update count
+bump :: Stats a -> Stats a
+bump s = s { total = total s + 1 }
+
+-- | Update data structure with a single alignment
+class1 :: Insertable x => Stats x -> Bam1 -> Stats x
+class1 c0 b
+  | isUnmapped b = c0
+  | isOpposite b =
+      (if firstUpstream b then add_innie else add_outie) b c0
+  | otherwise =
+      (if firstUpstream b then add_rightie else add_leftie) b c0
+
+add_innie, add_outie, add_rightie, add_leftie :: Insertable x => Bam1 -> Stats x -> Stats x
+add_innie b c = c { innies = insert b (innies c) }
+add_outie b c = c { outies = insert b (outies c) }
+add_rightie b c = c { righties = insert b (righties c) }
+add_leftie b c = c { lefties = insert b (lefties c) }
+
+isUnmapped, isOpposite, firstUpstream, isBefore :: Bam1 -> Bool
+isUnmapped = isNothing . insertSize
+-- isUnmap b || isMateUnmap b || mateTargetID b /= targetID b -- not the same
+-- insertSize is Nothing if negative (downstream read)
+isOpposite b = isReverse b /= isMateReverse b 
+firstUpstream b = isReverse b /= isBefore b
+isBefore b = position b < matePosition b
+
+-- --------------------------------------------------
+-- ClassStats for calculating mean, stdev, etc etc
+-- --------------------------------------------------
+
+data ClassStats = CS { ccount :: !Int, xsum, x2sum, x3sum, x4sum :: !Double } 
+                deriving Show
+
+data Statistics = S { ncount, mean, stdev, skew, kurt :: Double }
+
+mkstats :: ClassStats -> Statistics
+mkstats cs = S n m s sk kt
+  where
+    n = fromIntegral (ccount cs)
+    x = xsum cs
+    x2 = x2sum cs
+    x3 = x3sum cs
+    x4 = x4sum cs
+    m = x/n
+    m2 = m*m
+    m3 = m*m2
+    s = sqrt ((x2-m2*n)/(n-1))
+    sk = (x3 - 3*m*x2 + 2*m3*n)/(s*s*s*n)
+    kt = (x4 - 4*m*x3 + 6*m2*x2 - 4*m3*x + n*m*m3)/(s*s*s*s*n) - 3
+
+statistics :: Stats ClassStats -> Stats Statistics
+statistics (Class t i o l r) = (Class t (mkstats i) (mkstats o) (mkstats l) (mkstats r)) -- todo: Functor instance?
+
+instance Insertable ClassStats where
+  insert b (CS c s s2 s3 s4) = CS (c+1) (s+d) (s2+d^(2::Int)) (s3+d^(3::Int)) (s4+d^(4::Int)) 
+    where d = maybe (error ("no insert size?\n"++show b)) fromIntegral $ insertSize b
+  dispheader _ = ["         count","   prop","   mean","  stdev","   skew","   kurt"]
+  cdef = CS 0 0 0 0 0
+  disp1 tot cs = printf "%14d" (ccount cs) 
+               : printf "%5.2f%%" (200*fromIntegral (ccount cs)/fromIntegral tot::Double)
+               : map (printf "%7.1f") [mean s, stdev s, skew s, kurt s]
+    where 
+      s = mkstats cs  
+
+-- --------------------------------------------------
+-- Hist for accumulating histograms of sizes
+-- --------------------------------------------------
+
+data Hist = H { hcount :: !Int, buckets :: ![(Int,Int)] } deriving Show
+
+instance Insertable Hist where
+  insert b h = H (hcount h + 1) (go (fromIntegral $ fromJust $ insertSize b) (buckets h))
+    where go x ((b1,v):bs@(_:_)) =
+            if x > b1 then (b1,v):go x bs else inc (b1,v):bs
+          go x [(b1,v)] = [inc (b1,v)]
+          inc (x,y) = let y' = y+1 in y' `seq` (x,y')
+  disp1 tot h = printf "%14d" (hcount h)
+              : printf "%5.2f%%" (200*fromIntegral (hcount h)/fromIntegral tot::Double)
+              : map (printf "%7d" . snd) (buckets h)
+  dispheader h = "         count":"  prop":(map (printf "%7d" . fst) $ init $ buckets h)++["    >"]
+  cdef = undefined
+
+-- --------------------------------------------------
+-- Just "Collect" all the Bams in different classes
+-- --------------------------------------------------
+
+data Collect = Bams { bams :: [Bam1] }
+
+instance Insertable Collect where
+  insert b (Bams bs) = Bams (b:bs)
+  disp1 _ (Bams bs) = [unlines $ ("":map show bs)]
+  dispheader = const ["foo"]
+  cdef = Bams []