diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Marco Zocca (c) 2021
+
+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 Marco Zocca 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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,5 @@
+# rp-tree
+
+![scatterplot](/r/scatter.png "scatterplot")
+
+Random projection trees
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/app/Main.hs b/app/Main.hs
new file mode 100644
--- /dev/null
+++ b/app/Main.hs
@@ -0,0 +1,123 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE LambdaCase #-}
+{-# options_ghc -Wno-unused-imports #-}
+module Main where
+
+import Data.Foldable (fold)
+-- conduit
+import qualified Data.Conduit as C (ConduitT, runConduit, yield, await)
+import Data.Conduit ((.|))
+import qualified Data.Conduit.Combinators as C (map, mapM, scanl, scanlM, last, print)
+import qualified Data.Conduit.List as C (chunksOf, unfold, unfoldM)
+-- containers
+import qualified Data.IntMap as IM (IntMap, fromList, insert, lookup, map, mapWithKey, traverseWithKey, foldlWithKey, foldrWithKey)
+-- exceptions
+import Control.Monad.Catch (MonadThrow(..))
+-- splitmix-distributions
+import System.Random.SplitMix.Distributions (Gen, GenT, sample, sampleT, bernoulli, normal)
+-- transformers
+import Control.Monad.Trans.State.Lazy (State, get, put, evalState)
+-- vector
+import qualified Data.Vector as V (Vector, toList, fromList, replicate, zip)
+
+import Control.Monad (replicateM)
+import Data.RPTree (knn, candidates, Inner(..), RPTree, RPForest, leaves, SVector, DVector, fromListDv, dense, writeCsv, forest, dataSource)
+
+main :: IO ()
+main = do
+  let
+    n = 10000
+  -- renderTree0 n
+  -- renderTree1 n
+  undefined -- FIXME
+
+
+-- renderTree0 :: Int -> IO ()
+renderTree0 tt = do
+  let csvrows = V.toList $ fold $ flip evalState A $ traverse labeledV tt -- (tree0 n)
+  writeCsv "r/scatter_data.csv" csvrows
+
+-- renderTree1 :: Int -> IO ()
+renderTree1 tt = do
+  let
+    -- csvrows :: [(DVector Double, Pal5)]
+    csvrows = fold $ flip evalState A $ traverse labeledV tt -- (tree1 n)
+  writeCsv "r/scatter_data_rt2.csv" $ V.toList csvrows
+
+labeled :: (Enum i) =>
+           [a] -> State i [(a, i)]
+labeled xs = do
+  i <- get
+  put (succ i)
+  let n = length xs
+  pure $ zip xs (replicate n i)
+
+
+labeledV :: Enum i => V.Vector a -> State i (V.Vector (a, i))
+labeledV xs = do
+  i <- get
+  put (succ i)
+  let n = length xs
+  pure $ V.zip xs (V.replicate n i)
+
+data Pal5 = A | B | C | D | E deriving (Eq, Show)
+instance Enum Pal5 where
+  toEnum = \case
+    0 -> A
+    1 -> B
+    2 -> C
+    3 -> D
+    4 -> E
+    x -> toEnum (x `mod` 5)
+  fromEnum = \case
+    A -> 0
+    B -> 1
+    C -> 2
+    D -> 3
+    E -> 4
+
+-- tree0 :: Int -> RPTree Double (V.Vector (DVector Double))
+-- tree0 n = evalGen 1234 $ tree 10 1.0 2 (dataset n)
+
+-- tree1 :: Int -> RT SVector Double (V.Vector (DVector Double))
+-- tree1 n = evalGen 1234 $ treeRT 10 20 1.0 2 (dataset n)
+
+dataset :: Int -> V.Vector (DVector Double)
+dataset n = V.fromList $ sample 1234 $ replicateM n (dense 2 $ normal 0 1)
+
+
+-- treeC0 :: MonadThrow m =>
+--           Int -> GenT m (RPTree Double (V.Vector (DVector Double)))
+-- treeC0 n = treeSink 1234 10 20 100 1.0 2 (srcC n)
+
+{-
+λ> nn0 10000 (fromListDv [0,0])
+[(0.13092191004810114,DV [-8.771274989760332e-2,9.71957819868927e-2]),(0.14722273682679538,DV [-4.767722969780902e-2,0.13928896584839093]),(0.1626065099556818,DV [-4.57842765697381e-2,0.15602780873598454]),(0.22082909577433263,DV [-3.62336905451185e-2,0.21783619811681887]),(0.22085935710897311,DV [0.21196201255823421,-6.2056110535964756e-2]),(0.2636139991233282,DV [-0.24290511334764195,0.10241799862994452]),(0.3869415454995779,DV [-0.3658837577279577,0.12590804368455188]),(0.3951528583078011,DV [-0.3543713488257354,0.1748334308999686]),(0.6174219338196472,DV [-0.4952807707701239,0.3686553979897009]),(0.6968774335522048,DV [-0.6408548616154526,0.2737575638007956])]
+-}
+nn0 :: (Inner SVector v, Inner DVector v) =>
+       Int -> v Double -> V.Vector (Double, DVector Double)
+nn0 n q = case ttsm of
+  Just tts -> knn metricL2 10 tts q -- FIXME voting search size ?!
+  -- Nothing -> mempty
+  where
+    ttsm = sampleT 1234 $ forestC0 n
+
+cs0 n q = case sampleT 1234 $ forestC0 n of
+  Just tts -> (`candidates` q) <$> tts
+
+forestC0 :: MonadThrow m =>
+            Int
+         -> GenT
+            m
+            (IM.IntMap (RPTree Double (V.Vector (DVector Double))))
+forestC0 n = forest 1234 10 20 10 100 1.0 2 (srcC n)
+
+srcC :: Monad m => Int -> C.ConduitT i (DVector Double) (GenT m) ()
+srcC n = dataSource n normal2
+
+normal2 :: (Monad m) => GenT m (DVector Double)
+normal2 = do
+  b <- bernoulli 0.5
+  if b
+    then dense 2 $ normal 0 0.5
+    else dense 2 $ normal 2 0.5
diff --git a/rp-tree.cabal b/rp-tree.cabal
new file mode 100644
--- /dev/null
+++ b/rp-tree.cabal
@@ -0,0 +1,77 @@
+name:                rp-tree
+version:             0.1.0.0
+synopsis:            Random projection trees
+description:         Random projection trees for approximate nearest neighbor search in high-dimensional vector spaces
+homepage:            https://github.com/ocramz/rp-tree
+license:             BSD3
+license-file:        LICENSE
+author:              Marco Zocca
+maintainer:          ocramz
+copyright:           2021 Marco Zocca
+category:            Data Mining, Data Structures, Machine Learning, Data
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+tested-with:         GHC == 8.10.4
+
+library
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  hs-source-dirs:      src
+  exposed-modules:     Data.RPTree
+
+  other-modules:       Data.RPTree.Internal
+                       Data.RPTree.Gen
+                       Data.RPTree.Draw
+                       Data.RPTree.Conduit
+  build-depends:       base >= 4.7 && < 5
+                     , boxes
+                     , bytestring
+                     , conduit
+                     , containers >= 0.6.2.1
+                     , deepseq >= 1.4.4.0
+                     , exceptions
+                     , microlens
+                     , microlens-th
+                     , mtl
+                     -- , psqueues
+                     , serialise
+                     , splitmix-distributions >= 0.8
+                     , transformers
+                     -- , ulid
+                     , vector >= 0.12.1.2
+                     , vector-algorithms
+                       -- debug
+                     , hspec
+
+test-suite spec
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  build-depends:       base
+                     , rp-tree
+                     , hspec
+                     , QuickCheck
+                     , splitmix-distributions
+
+executable rp-tree
+  default-language:    Haskell2010
+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+  hs-source-dirs:      app
+  main-is:             Main.hs
+  build-depends:       base
+                     , conduit
+                     , containers
+                     , exceptions
+                     , rp-tree
+                     , splitmix-distributions
+                     , transformers
+                     , vector
+
+source-repository head
+  type:     git
+  location: https://github.com/ocramz/rp-tree
+
+
diff --git a/src/Data/RPTree.hs b/src/Data/RPTree.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/RPTree.hs
@@ -0,0 +1,359 @@
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# language DeriveGeneric #-}
+{-# language LambdaCase #-}
+{-# language GeneralizedNewtypeDeriving #-}
+-- {-# language MultiParamTypeClasses #-}
+{-# LANGUAGE MultiWayIf #-}
+{-# options_ghc -Wno-unused-imports #-}
+{-# options_ghc -Wno-unused-top-binds #-}
+
+{-|
+Random projection trees for approximate nearest neighbor search in high-dimensional vector spaces
+-}
+module Data.RPTree (
+  -- * Construction
+  forest
+  -- * Query
+  , knn
+  -- , nearest
+  -- * Validation
+  , recall
+  -- * Access
+  , levels, points, leaves, candidates
+  -- * Types
+  -- ** RPTree
+  , RPTree, RPForest
+  -- -- *** internal
+  -- , RPT
+  -- -- ** RT
+  -- , RT
+  -- *
+  , SVector, fromListSv
+  , DVector, fromListDv
+  -- * inner product
+  , Inner(..), Scale(..)
+  --   -- ** helpers for implementing Inner instances
+  --   -- *** inner product
+  -- , innerSS, innerSD
+  --   -- *** L2 distance
+  -- , metricSSL2, metricSDL2
+  -- * Conduit
+  , dataSource
+  -- * Random generation
+  -- ** vector
+  , sparse, dense
+  -- * Rendering
+  , draw
+  -- * CSV
+  , writeCsv
+  ) where
+
+import Control.Monad (replicateM)
+
+import Control.Monad.IO.Class (MonadIO(..))
+import Data.Foldable (Foldable(..), maximumBy, minimumBy)
+import Data.Functor.Identity (Identity(..))
+import Data.List (partition, sortBy)
+import Data.Monoid (Sum(..))
+import Data.Ord (comparing)
+import GHC.Generics (Generic)
+import GHC.Word (Word64)
+
+-- containers
+import Data.Sequence (Seq, (|>))
+import qualified Data.Map as M (Map, fromList, toList, foldrWithKey, insert, insertWith)
+import qualified Data.Set as S (Set, fromList, intersection, insert)
+-- deepseq
+import Control.DeepSeq (NFData(..))
+-- mtl
+import Control.Monad.State (MonadState(..), modify)
+-- -- psqueues
+-- import qualified Data.IntPSQ as PQ (IntPSQ, insert, fromList, findMin, minView)
+-- transformers
+import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState)
+import Control.Monad.Trans.Class (MonadTrans(..))
+-- vector
+import qualified Data.Vector as V (Vector, replicateM, fromList)
+import qualified Data.Vector.Generic as VG (Vector(..), unfoldrM, length, replicateM, (!), map, freeze, thaw, take, drop, unzip)
+import qualified Data.Vector.Unboxed as VU (Vector, Unbox, fromList)
+import qualified Data.Vector.Storable as VS (Vector)
+-- vector-algorithms
+import qualified Data.Vector.Algorithms.Merge as V (sortBy)
+
+import Data.RPTree.Conduit (forest, dataSource)
+import Data.RPTree.Gen (sparse, dense)
+import Data.RPTree.Internal (RPTree(..), RPForest, RPT(..), levels, points, leaves, RT(..), Inner(..), Scale(..), (/.), innerSD, innerSS, metricSSL2, metricSDL2, SVector(..), fromListSv, DVector(..), fromListDv, partitionAtMedian, Margin, getMargin, sortByVG)
+
+import Data.RPTree.Draw (draw, writeCsv)
+
+
+-- | k nearest neighbors
+knn :: (Ord p, Inner SVector v, VU.Unbox d, Real d) =>
+       (v2 -> v d -> p) -- ^ distance function
+    -> Int -- ^ k neighbors
+    -> RPForest d (V.Vector v2) -- ^ random projection forest
+    -> v d -- ^ query point
+    -> V.Vector (p, v2) -- ^ ordered in increasing distance order
+knn distf k tts q = sortByVG fst cs
+  where
+    cs = VG.map (\x -> (x `distf` q, x)) $ VG.take k $ fold $ (`candidates` q) <$> tts
+
+
+-- | average recall-at-k, computed over a set of trees
+recall :: (Inner u v, Inner SVector v, VU.Unbox a, Ord a,
+            Ord (u a), Floating a) =>
+          RPForest a (V.Vector (u a))
+       -> Int -- ^ k : number of nearest neighbors to consider
+       -> v a -- ^ query point
+       -> Double
+recall tt k q = sum rs / fromIntegral n
+  where
+    rs = fmap (\t -> recall1 t k q) tt
+    n = length tt
+
+recall1 :: (Inner SVector v, Inner u v, VU.Unbox a, Ord a, Ord (u a), Floating a) =>
+          RPTree a (V.Vector (u a))
+       -> Int -- ^ k : number of nearest neighbors to consider
+       -> v a  -- ^ query point
+       -> Double
+recall1 = recallWith metricL2
+
+recallWith :: (Fractional a1, Inner SVector v, Ord d, VU.Unbox d,
+                Num d, Ord a3, Ord a2) =>
+              (a2 -> v d -> a3) -> RPTree d (V.Vector a2) -> Int -> v d -> a1
+recallWith distf tt k q = fromIntegral (length aintk) / fromIntegral k
+  where
+    xs = points tt
+    dists = sortBy (comparing snd) $ toList $ fmap (\x -> (x, x `distf` q)) xs
+    kk = S.fromList $ map fst $ take k dists
+    aa = set $ candidates tt q
+    aintk = aa `S.intersection` kk
+
+set :: (Foldable t, Ord a) => t a -> S.Set a
+set = foldl (flip S.insert) mempty
+
+
+
+-- | Retrieve points nearest to the query
+--
+-- in case of a narrow margin, collect both branches of the tree
+candidates :: (Inner SVector v, VU.Unbox d, Ord d, Num d, Semigroup xs) =>
+              RPTree d xs
+           -> v d -- ^ query point
+           -> xs
+candidates (RPTree rvs tt) x = go 0 tt
+  where
+    go _     (Tip xs)                     = xs
+    go ixLev (Bin thr margin ltree rtree) = do
+      let
+        (mglo, mghi) = getMargin margin
+        r = rvs VG.! ixLev
+        proj = r `inner` x
+        i' = succ ixLev
+      if | proj < thr &&
+           mglo > mghi -> go i' ltree <> go i' rtree
+         | proj < thr  -> go i' ltree
+         | proj > thr &&
+           mglo < mghi -> go i' ltree <> go i' rtree
+         | otherwise   -> go i' rtree
+
+
+
+
+
+
+-- pqSeq :: Ord a => PQ.IntPSQ a b -> Seq (a, b)
+-- pqSeq pqq = go pqq mempty
+--   where
+--     go pq acc = case PQ.minView pq of
+--       Nothing -> acc
+--       Just (_, p, v, rest) -> go rest $ acc |> (p, v)
+
+
+-- newtype Counts a = Counts {
+--   unCounts :: M.Map a (Sum Int) } deriving (Eq, Show, Semigroup, Monoid)
+-- keepCounts :: Int -- ^ keep entry iff counts are larger than this value
+--            -> Counts a
+--            -> [(a, Int)]
+-- keepCounts thr cs = M.foldrWithKey insf mempty c
+--   where
+--     insf k v acc
+--       | v >= thr = (k, v) : acc
+--       | otherwise = acc
+--     c = getSum `fmap` unCounts cs
+-- counts :: (Foldable t, Ord a) => t a -> Counts a
+-- counts = foldl count mempty
+-- count :: Ord a => Counts a -> a -> Counts a
+-- count (Counts mm) x = Counts $ M.insertWith mappend x (Sum 1) mm
+
+
+-- forest :: Inner SVector v =>
+--           Int -- ^ # of trees
+--        -> Int -- ^ maximum tree height
+--        -> Double -- ^ nonzero density of sparse projection vectors
+--        -> Int -- ^ dimension of projection vectors
+--        -> V.Vector (v Double) -- ^ dataset
+--        -> Gen [RPTree Double (V.Vector (v Double))]
+-- forest nt maxDepth pnz dim xss =
+--   replicateM nt (tree maxDepth pnz dim xss)
+
+-- -- | Build a random projection tree
+-- --
+-- -- Optimization: instead of sampling one projection vector per branch, we sample one per tree level (as suggested in https://www.cs.helsinki.fi/u/ttonteri/pub/bigdata2016.pdf )
+-- tree :: (Inner SVector v) =>
+--          Int -- ^ maximum tree height
+--       -> Double -- ^ nonzero density of sparse projection vectors
+--       -> Int -- ^ dimension of projection vectors
+--       -> V.Vector (v Double) -- ^ dataset
+--       -> Gen (RPTree Double (V.Vector (v Double)))
+-- tree maxDepth pnz dim xss = do
+--   -- sample all projection vectors
+--   rvs <- V.replicateM maxDepth (sparse pnz dim stdNormal)
+--   let
+--     loop ixLev xs = do
+--       if ixLev >= maxDepth || length xs <= 100
+--         then
+--           pure $ Tip xs
+--         else
+--         do
+--           let
+--             r = rvs VG.! ixLev
+--             (thr, margin, ll, rr) = partitionAtMedian r xs
+--           treel <- loop (ixLev + 1) ll
+--           treer <- loop (ixLev + 1) rr
+--           pure $ Bin thr margin treel treer
+--   rpt <- loop 0 xss
+--   pure $ RPTree rvs rpt
+
+
+
+
+
+-- -- | Partition at median inner product
+-- treeRT :: (Monad m, Inner SVector v) =>
+--            Int
+--         -> Int
+--         -> Double
+--         -> Int
+--         -> V.Vector (v Double)
+--         -> GenT m (RT SVector Double (V.Vector (v Double)))
+-- treeRT maxDepth minLeaf pnz dim xss = loop 0 xss
+--   where
+--     loop ixLev xs = do
+--       if ixLev >= maxDepth || length xs <= minLeaf
+--         then
+--           pure $ RTip xs
+--         else
+--         do
+--           r <- sparse pnz dim stdNormal
+--           let
+--             (_, mrg, ll, rr) = partitionAtMedian r xs
+--           treel <- loop (ixLev + 1) ll
+--           treer <- loop (ixLev + 1) rr
+--           pure $ RBin r mrg treel treer
+
+
+
+
+
+
+
+-- -- | Like 'tree' but here we partition at the median of the inner product values instead
+-- tree' :: (Inner SVector v) =>
+--          Int
+--       -> Double
+--       -> Int
+--       -> V.Vector (v Double)
+--       -> Gen (RPTree Double (V.Vector (v Double)))
+-- tree' maxDepth pnz dim xss = do
+--   -- sample all projection vectors
+--   rvs <- V.replicateM maxDepth (sparse pnz dim stdNormal)
+--   let
+--     loop ixLev xs =
+--       if ixLev >= maxDepth || length xs <= 100
+--         then Tip xs
+--         else
+--           let
+--             r = rvs VG.! ixLev
+--             (thr, margin, ll, rr) = partitionAtMedian r xs
+--             tl = loop (ixLev + 1) ll
+--             tr = loop (ixLev + 1) rr
+--           in Bin thr margin tl tr
+--   let rpt = loop 0 xss
+--   pure $ RPTree rvs rpt
+
+
+-- -- | Partition uniformly at random between inner product extreme values
+-- treeRT :: (Monad m, Inner SVector v) =>
+--           Int -- ^ max tree depth
+--        -> Int -- ^ min leaf size
+--        -> Double -- ^ nonzero density
+--        -> Int -- ^ embedding dimension
+--        -> V.Vector (v Double) -- ^ data
+--        -> GenT m (RT SVector Double (V.Vector (v Double)))
+-- treeRT maxDepth minLeaf pnz dim xss = loop 0 xss
+--   where
+--     loop ixLev xs = do
+--       if ixLev >= maxDepth || length xs <= minLeaf
+--         then
+--           pure $ RTip xs
+--         else
+--         do
+--           -- sample projection vector
+--           r <- sparse pnz dim stdNormal
+--           let
+--             -- project the dataset
+--             projs = map (\x -> (x, r `inner` x)) xs
+--             hi = snd $ maximumBy (comparing snd) projs
+--             lo = snd $ minimumBy (comparing snd) projs
+--           -- sample a threshold
+--           thr <- uniformR lo hi
+--           let
+--             (ll, rr) = partition (\xp -> snd xp < thr) projs
+--           treel <- loop (ixLev + 1) (map fst ll)
+--           treer <- loop (ixLev + 1) (map fst rr)
+--           pure $ RBin r treel treer
+
+
+-- -- | Partition wrt a plane _|_ to the segment connecting two points sampled at random
+-- --
+-- -- (like annoy@@)
+-- treeRT2 :: (Monad m, Ord d, Fractional d, Inner v v, VU.Unbox d, Num d) =>
+--            Int
+--         -> Int
+--         -> [v d]
+--         -> GenT m (RT v d [v d])
+-- treeRT2 maxd minl xss = loop 0 xss
+--   where
+--     loop ixLev xs = do
+--       if ixLev >= maxd || length xs <= minl
+--         then
+--           pure $ RTip xs
+--         else
+--         do
+--           x12 <- sampleWOR 2 xs
+--           let
+--             (x1:x2:_) = x12
+--             r = x1 ^-^ x2
+--             (ll, rr) = partition (\x -> (r `inner` (x ^-^ x1) < 0)) xs
+--           treel <- loop (ixLev + 1) ll
+--           treer <- loop (ixLev + 1) rr
+--           pure $ RBin r treel treer
+
+
+
+
+
+
+
+
+
+
+-- ulid :: MonadIO m => a -> m (ULID a)
+-- ulid x = ULID <$> pure x <*> liftIO UU.getULID
+-- data ULID a = ULID { uData :: a , uULID :: UU.ULID } deriving (Eq, Show)
+-- instance (Eq a) => Ord (ULID a) where
+--   ULID _ u1 <= ULID _ u2 = u1 <= u2
diff --git a/src/Data/RPTree/Conduit.hs b/src/Data/RPTree/Conduit.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/RPTree/Conduit.hs
@@ -0,0 +1,215 @@
+{-# language DeriveDataTypeable #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# options_ghc -Wno-unused-imports #-}
+{-# options_ghc -Wno-unused-top-binds #-}
+module Data.RPTree.Conduit (
+  forest
+  -- ** helpers
+  , dataSource
+  ) where
+
+import Control.Monad (replicateM)
+import GHC.Word (Word64)
+
+-- conduit
+import qualified Data.Conduit as C (ConduitT, runConduit, yield, await)
+import Data.Conduit ((.|))
+import qualified Data.Conduit.Combinators as C (map, mapM, scanl, scanlM, last, print)
+import qualified Data.Conduit.List as C (chunksOf, unfold, unfoldM)
+-- containers
+import qualified Data.IntMap as IM (IntMap, fromList, insert, lookup, map, mapWithKey, traverseWithKey, foldlWithKey, foldrWithKey, intersectionWith)
+-- exceptions
+import Control.Monad.Catch (MonadThrow(..))
+-- mtl
+import Control.Monad.State (MonadState(..), modify)
+-- splitmix-distributions
+import System.Random.SplitMix.Distributions (Gen, sample, GenT, sampleT, normal, stdNormal, stdUniform, exponential, bernoulli, uniformR)
+-- transformers
+import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState)
+import Control.Monad.Trans.Class (MonadTrans(..))
+-- vector
+import qualified Data.Vector as V (Vector, replicateM, fromList)
+import qualified Data.Vector.Generic as VG (Vector(..), unfoldrM, length, replicateM, (!), map, freeze, thaw, take, drop, unzip)
+import qualified Data.Vector.Unboxed as VU (Vector, Unbox, fromList)
+import qualified Data.Vector.Storable as VS (Vector)
+
+import Data.RPTree.Gen (sparse, dense)
+import Data.RPTree.Internal (RPTree(..), RPForest, RPT(..), levels, points, Inner(..), innerSD, innerSS, metricSSL2, metricSDL2, SVector(..), fromListSv, DVector(..), fromListDv, partitionAtMedian, RPTError(..))
+
+
+
+
+-- | Source of random data points
+dataSource :: (Monad m) =>
+              Int -- ^ number of vectors to generate
+           -> GenT m a -- ^ random generator for the vector components
+           -> C.ConduitT i a (GenT m) ()
+dataSource n gg = flip C.unfoldM 0 $ \i -> do
+  if i == n
+    then pure Nothing
+    else do
+      x <- gg
+      pure $ Just (x, i + 1)
+
+-- | Populate a tree from a data stream
+--
+-- Assumptions on the data source:
+--
+-- * non-empty : contains at least one value
+--
+-- * stationary : each chunk is representative of the whole dataset
+--
+-- * bounded : we wait until the end of the stream to produce a result
+--
+-- Throws 'EmptyResult' if the conduit is empty
+tree :: (MonadThrow m, Inner SVector v) =>
+            Word64 -- ^ random seed
+         -> Int -- ^ max tree depth
+         -> Int -- ^ min leaf size
+         -> Int -- ^ data chunk size
+         -> Double -- ^ nonzero density of projection vectors
+         -> Int -- ^ dimension of projection vectors
+         -> C.ConduitT () (v Double) m () -- ^ data source
+         -> m (RPTree Double (V.Vector (v Double)))
+tree seed maxDepth minLeaf n pnz dim src = do
+  let
+    rvs = sample seed $ V.replicateM maxDepth (sparse pnz dim stdNormal)
+  tm <- C.runConduit $ src .|
+                       insertC maxDepth minLeaf n rvs .|
+                       C.last
+  case tm of
+    Just t -> pure $ RPTree rvs t
+    _ -> throwM $ EmptyResult "treeSink"
+
+-- | Incrementally build a tree
+insertC :: (Monad m, Inner u v, Ord d, VU.Unbox d, Fractional d) =>
+           Int -- ^ max tree depth
+        -> Int -- ^ min leaf size
+        -> Int -- ^ data chunk size
+        -> V.Vector (u d) -- ^ random projection vectors
+        -> C.ConduitT (v d) (RPT d (V.Vector (v d))) m ()
+insertC maxDepth minLeaf n rvs = chunked n z (insert maxDepth minLeaf rvs)
+  where
+    z = Tip mempty
+
+
+
+-- | Populate a forest from a data stream
+--
+-- Assumptions on the data source:
+--
+-- * non-empty : contains at least one value
+--
+-- * stationary : each chunk is representative of the whole dataset
+--
+-- * bounded : we wait until the end of the stream to produce a result
+--
+-- Throws 'EmptyResult' if the conduit is empty
+forest :: (MonadThrow m, Inner SVector v) =>
+                 Word64 -- ^ random seed
+              -> Int -- ^ max tree depth
+              -> Int -- ^ min leaf size
+              -> Int -- ^ number of trees
+              -> Int -- ^ data chunk size
+              -> Double -- ^ nonzero density of projection vectors
+              -> Int -- ^ dimension of projection vectors
+              -> C.ConduitT () (v Double) m () -- ^ data source
+              -> m (RPForest Double (V.Vector (v Double)))
+forest seed maxd minl ntrees chunksize pnz dim src = do
+  let
+    rvss = sample seed $ do
+      rvs <- replicateM ntrees $ V.replicateM maxd (sparse pnz dim stdNormal)
+      pure $ IM.fromList $ zip [0 .. ] rvs
+  tm <- C.runConduit $ src .|
+                       insertMultiC maxd minl chunksize rvss .|
+                       C.last
+  case tm of
+    Just ts -> do
+      let
+        res = IM.intersectionWith RPTree rvss ts
+      pure res
+    _ -> throwM $ EmptyResult "forestSink"
+
+
+
+
+
+insertMultiC :: (Monad m, Ord d, Inner u v, VU.Unbox d, Fractional d, VG.Vector v1 (u d)) =>
+                Int  -- ^ max tree depth
+             -> Int -- ^ min leaf size
+             -> Int -- ^ chunk size
+             -> IM.IntMap (v1 (u d)) -- one entry per tree
+             -> C.ConduitT
+                (v d)
+                (IM.IntMap (RPT d (V.Vector (v d))))
+                m
+                ()
+insertMultiC maxd minl n rvss = chunked n im0 (insertMulti maxd minl rvss)
+  where
+    im0 = IM.map (const z) rvss
+    z = Tip mempty
+
+
+insertMulti :: (Ord d, Inner u v, VU.Unbox d, Fractional d, VG.Vector v1 (u d)) =>
+               Int
+            -> Int
+            -> IM.IntMap (v1 (u d)) -- ^ projection vectors
+            -> IM.IntMap (RPT d (V.Vector (v d))) -- ^ accumulator of subtrees
+            -> V.Vector (v d) -- ^ data chunk
+            -> IM.IntMap (RPT d (V.Vector (v d)))
+insertMulti maxd minl rvss tacc xs =
+  flip IM.mapWithKey tacc $ \i t -> case IM.lookup i rvss of
+                                      Just rvs -> insert maxd minl rvs t xs
+                                      _        -> t
+
+insert :: (VG.Vector v1 (u d), Ord d, Inner u v, VU.Unbox d, Fractional d) =>
+          Int -- ^ max tree depth
+       -> Int -- ^ min leaf size
+       -> v1 (u d) -- ^ projection vectors
+       -> RPT d (V.Vector (v d)) -- ^ accumulator
+       -> V.Vector (v d) -- ^ data chunk
+       -> RPT d (V.Vector (v d))
+insert maxDepth minLeaf rvs = loop 0
+  where
+    z = Tip mempty
+    loop ixLev tt xs =
+      let
+        r = rvs VG.! ixLev
+      in
+        case tt of
+
+          b@(Bin thr0 margin0 tl0 tr0) ->
+            if ixLev >= maxDepth || length xs <= minLeaf
+              then b -- return current subtree
+              else
+                let
+                  (thr, margin, ll, rr) =
+                    partitionAtMedian r xs
+                  margin' = margin0 <> margin
+                  thr' = (thr0 + thr) / 2
+                  tl = loop (ixLev + 1) tl0 ll
+                  tr = loop (ixLev + 1) tr0 rr
+                in Bin thr' margin' tl tr
+
+          Tip xs0 -> do
+            let xs' = xs <> xs0
+            if ixLev >= maxDepth || length xs <= minLeaf
+              then Tip xs' -- concat data in leaf
+              else
+                let
+                  (thr, margin, ll, rr) = partitionAtMedian r xs'
+                  tl = loop (ixLev + 1) z ll
+                  tr = loop (ixLev + 1) z rr
+                in Bin thr margin tl tr
+
+
+chunked :: (Monad m) =>
+           Int -- ^ chunk size
+        -> t -- ^ initial tree
+        -> (t -> V.Vector a -> t)
+        -> C.ConduitT a t m ()
+chunked n z f = do C.chunksOf n .|
+                     C.map V.fromList .|
+                     C.scanl f z -- .|
+
+
diff --git a/src/Data/RPTree/Draw.hs b/src/Data/RPTree/Draw.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/RPTree/Draw.hs
@@ -0,0 +1,86 @@
+{-# language LambdaCase #-}
+{-# options_ghc -Wno-unused-imports #-}
+module Data.RPTree.Draw where
+
+import Data.List (intercalate)
+import Text.Printf (PrintfArg, printf)
+
+-- boxes
+import qualified Text.PrettyPrint.Boxes as B (Box, render, emptyBox, vcat, hcat, text, top, bottom, center1)
+-- bytestring
+import qualified Data.ByteString.Lazy    as LBS (ByteString, writeFile)
+import qualified Data.ByteString.Builder as BSB (Builder, toLazyByteString, string7, charUtf8)
+-- -- mtl
+-- import Control.Monad.State (MonadState(..))
+-- vector
+import qualified Data.Vector as V (Vector, replicateM)
+import qualified Data.Vector.Generic as VG (Vector(..), map, sum, unfoldr, unfoldrM, length, replicateM, (!))
+import qualified Data.Vector.Unboxed as VU (Unbox)
+
+import Data.RPTree.Internal (RPTree(..), RPT(..), DVector, toListDv)
+
+
+
+
+-- | Encode dataset as CSV and save into file
+writeCsv :: (Show a, Show b, VU.Unbox a) =>
+            FilePath
+         -> [(DVector a, b)] -- ^ data point, label
+         -> IO ()
+writeCsv fp ds = LBS.writeFile fp $ BSB.toLazyByteString $ toCsv ds
+
+toCsvRow :: (Show a, Show b, VU.Unbox a) =>
+            DVector a
+         -> b
+         -> BSB.Builder
+toCsvRow dv i = BSB.string7 $ intercalate "," [show x, show y, show i]
+  where
+    (x:y:_) = toListDv dv
+
+toCsv :: (Show a, Show b, VU.Unbox a) =>
+         [(DVector a, b)] -> BSB.Builder
+toCsv rs = mconcat [toCsvRow r i <> BSB.charUtf8 '\n' | (r, i) <- rs]
+
+-- | Render a tree to stdout
+--
+-- Useful for debugging
+--
+-- This should be called only for small trees, otherwise the printed result quickly overflows the screen and becomes hard to read.
+--
+-- NB : prints distance information rounded to two decimal digits
+draw :: (Show a, Boxed a, PrintfArg v) => RPTree v a -> IO ()
+draw = drawRPT . _rpTree
+
+drawRPT :: (Show a, Boxed a, PrintfArg v) => RPT v a -> IO ()
+drawRPT = putStrLn . toStringRPT
+
+toStringRPT :: (Show a, Boxed a, PrintfArg v) => RPT v a -> String
+toStringRPT = B.render . toBox
+
+toBox :: (Show a, Boxed a, PrintfArg v) => RPT v a -> B.Box
+toBox = \case
+  (Bin thr _ tl tr) ->
+    txt (node thr) `stack` (toBox tl `byside` toBox tr)
+  Tip xs -> boxed xs -- tipData xs -- txt $ show x
+  where
+    node x = printf "%5.2f" x -- (show x)
+
+class Boxed a where
+  boxed :: a -> B.Box
+instance (Show a) => Boxed [a] where
+  boxed = foldl (\bx x -> bx `stack` txt (show x)) $ B.emptyBox 0 0
+instance Boxed () where
+  boxed _ = txt "*"
+
+tipData :: (Show a, Foldable t) => t a -> B.Box
+tipData = foldl (\bx x -> bx `stack` txt (show x)) $ B.emptyBox 1 1
+
+txt :: String -> B.Box
+txt t = spc `byside` B.text t `byside` spc
+  where spc = B.emptyBox 1 1
+
+byside :: B.Box -> B.Box -> B.Box
+byside l r = B.hcat B.top [l, r]
+
+stack :: B.Box -> B.Box -> B.Box
+stack t b = B.vcat B.center1 [t, b]
diff --git a/src/Data/RPTree/Gen.hs b/src/Data/RPTree/Gen.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/RPTree/Gen.hs
@@ -0,0 +1,149 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# language LambdaCase #-}
+{-# language GeneralizedNewtypeDeriving #-}
+{-# options_ghc -Wno-unused-imports #-}
+module Data.RPTree.Gen where
+
+import Control.Monad (replicateM, foldM)
+
+-- containers
+import qualified Data.IntMap as IM (IntMap, insert, toList)
+-- mtl
+import Control.Monad.Trans.Class (MonadTrans(..))
+import Control.Monad.State (MonadState(..), modify)
+-- splitmix-distribitions
+import System.Random.SplitMix.Distributions (Gen, GenT, stdUniform, bernoulli)
+-- transformers
+import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState)
+-- vector
+
+
+import qualified Data.Vector.Generic as VG (Vector(..), unfoldrM, length, replicateM, (!))
+import qualified Data.Vector.Unboxed as VU (Vector, Unbox, fromList)
+
+
+import Data.RPTree.Internal (RPTree(..), RPT(..), SVector(..), fromListSv, DVector(..))
+
+
+-- | Sample without replacement with a single pass over the data
+--
+-- implements Algorithm L for reservoir sampling
+--
+-- Li, Kim-Hung (4 December 1994). "Reservoir-Sampling Algorithms of Time Complexity O(n(1+log(N/n)))". ACM Transactions on Mathematical Software. 20 (4): 481–493. doi:10.1145/198429.198435
+sampleWOR :: (Monad m, Foldable t) =>
+             Int -- ^ sample size
+          -> t a
+          -> GenT m [a]
+sampleWOR k xs = do
+  (_, res) <- flip runStateT z $ foldM insf 0 xs
+  pure $ map snd $ IM.toList (rsReservoir res)
+  where
+    z = RSPartial mempty
+    insf i x = do
+      st <- get
+      case st of
+        RSPartial acc -> do
+          w <- lift $ genW k
+          s <- lift $ genS w
+          let
+            acc' = IM.insert i x acc
+            ila = i + s + 1
+            st'
+              | i >= k = RSFull acc' ila w
+              | otherwise = RSPartial acc'
+          put st'
+          pure (succ i)
+        RSFull acc ila0 w0 -> do
+          case i `compare` ila0 of
+            EQ -> do
+              w <- lift $ genW k
+              s <- lift $ genS w0
+              let
+                ila = i + s + 1
+              acc' <- lift $ replaceInBuffer k acc x
+              let
+                w' = w0 * w
+              put (RSFull acc' ila w')
+              pure (succ i)
+            _ -> pure (succ i)
+
+data ResS a = RSPartial { rsReservoir :: IM.IntMap a }
+            | RSFull {
+                rsReservoir :: IM.IntMap a -- ^ reservoir
+                , rsfLookAh :: !Int -- ^ lookahead index
+                , rsfW :: !Double -- ^ W
+                } deriving (Eq, Show)
+
+genW :: (Monad m) => Int -> GenT m Double
+genW k = do
+  u <- stdUniform
+  pure $ exp (log u / fromIntegral k)
+
+genS :: (Monad m) => Double -> GenT m Int
+genS w = do
+  u <- stdUniform
+  pure $ floor (log u / log (1 - w))
+
+-- | Replaces a value at a random position within the buffer
+replaceInBuffer :: (Monad m) =>
+                   Int
+                -> IM.IntMap a
+                -> a
+                -> GenT m (IM.IntMap a)
+replaceInBuffer k imm y = do
+  u <- stdUniform
+  let ix = floor (fromIntegral k * u)
+  pure $ IM.insert ix y imm
+
+
+
+-- | Generate a sparse random vector with a given nonzero density and components sampled from the supplied random generator
+sparse :: (Monad m, VU.Unbox a) =>
+          Double -- ^ nonzero density
+       -> Int -- ^ vector dimension
+       -> GenT m a -- ^ random generator of vector components
+       -> GenT m (SVector a)
+sparse p sz rand = SV sz <$> sparseVG p sz rand
+
+-- | Generate a dense random vector with components sampled from the supplied random generator
+dense :: (Monad m, VG.Vector VU.Vector a) =>
+         Int -- ^ vector dimension
+      -> GenT m a -- ^ random generator of vector components
+      -> GenT m (DVector a)
+dense sz rand = DV <$> denseVG sz rand
+
+
+
+-- | Sample a dense random vector
+denseVG :: (VG.Vector v a, Monad m) =>
+           Int -- ^ vector dimension
+        -> m a
+        -> m (v a)
+denseVG sz rand = VG.unfoldrM mkf 0
+  where
+    mkf i
+      | i >= sz = pure Nothing
+      | otherwise = do
+          x <- rand
+          pure $ Just (x, succ i)
+
+-- | Sample a sparse random vector
+sparseVG :: (Monad m, VG.Vector v (Int, a)) =>
+            Double -- ^ nonzero density
+         -> Int  -- ^ vector dimension
+         -> GenT m a
+         -> GenT m (v (Int, a))
+sparseVG p sz rand = VG.unfoldrM mkf 0
+  where
+    mkf i
+      | i >= sz = pure Nothing
+      | otherwise = do
+          flag <- bernoulli p
+          if flag
+            then
+            do
+              x <- rand
+              pure $ Just ((i, x), succ i)
+            else
+              mkf (succ i)
diff --git a/src/Data/RPTree/Internal.hs b/src/Data/RPTree/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/RPTree/Internal.hs
@@ -0,0 +1,380 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE DeriveTraversable #-}
+{-# LANGUAGE DeriveFoldable #-}
+-- {-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# language DeriveGeneric #-}
+{-# language LambdaCase #-}
+{-# language MultiParamTypeClasses #-}
+{-# language GeneralizedNewtypeDeriving #-}
+{-# language TemplateHaskell #-}
+{-# options_ghc -Wno-unused-imports #-}
+module Data.RPTree.Internal where
+
+import Control.Exception (Exception(..))
+import Control.Monad.IO.Class (MonadIO(..))
+import Control.Monad.ST (runST)
+import Data.Function ((&))
+import Data.Foldable (fold, foldl')
+import Data.Functor.Identity (Identity(..))
+import Data.List (nub)
+import Data.Monoid (Sum(..))
+import Data.Ord (comparing)
+import Data.Semigroup (Min(..), Max(..))
+import Data.Typeable (Typeable)
+import GHC.Generics (Generic)
+
+-- containers
+import qualified Data.IntMap as IM (IntMap)
+-- deepseq
+import Control.DeepSeq (NFData(..))
+-- microlens
+import Lens.Micro (Traversal', (.~), (^..), folded)
+import Lens.Micro.TH (makeLensesFor, makeLensesWith, lensRules, generateSignatures)
+-- mtl
+import Control.Monad.State (MonadState(..), modify)
+-- serialise
+import Codec.Serialise (Serialise(..))
+-- transformers
+import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState)
+-- vector
+import qualified Data.Vector as V (Vector, replicateM, fromList)
+import qualified Data.Vector.Generic as VG (Vector(..), map, sum, unfoldr, unfoldrM, length, replicateM, (!), take, drop, unzip, freeze, thaw, foldl, foldr, toList, zipWith, last, head)
+import qualified Data.Vector.Unboxed as VU (Vector, Unbox, fromList, toList)
+import qualified Data.Vector.Storable as VS (Vector)
+-- vector-algorithms
+import qualified Data.Vector.Algorithms.Merge as V (sortBy)
+
+
+-- | Exceptions
+data RPTError =
+  EmptyResult String
+  deriving (Eq, Typeable)
+instance Show RPTError where
+  show = \case
+    EmptyResult str -> unwords [str, ": empty result"]
+instance Exception RPTError
+
+-- | Bounds around the cutting plane
+data Margin a = Margin {
+  cMarginLow :: Max a -- ^ lower bound on the cut point
+  , cMarginHigh :: Min a -- ^ upper bound
+                   } deriving (Eq, Show, Generic)
+instance (Serialise a) => Serialise (Margin a)
+getMargin :: Margin a -> (a, a)
+getMargin (Margin ml mh) = (getMax ml, getMin mh)
+instance (NFData a) => NFData (Margin a)
+-- | Used for updating in a streaming setting
+instance (Ord a) => Semigroup (Margin a) where
+  Margin lo1 hi1 <> Margin lo2 hi2 = Margin (lo1 <> lo2) (hi1 <> hi2)
+
+
+-- | Sparse vectors with unboxed components
+data SVector a = SV { svDim :: !Int, svVec :: VU.Vector (Int, a) } deriving (Eq, Ord, Generic)
+instance (VU.Unbox a, Serialise a) => Serialise (SVector a)
+instance (VU.Unbox a, Show a) => Show (SVector a) where
+  show (SV n vv) = unwords ["SV", show n, show (VU.toList vv)]
+instance NFData (SVector a)
+
+fromListSv :: VU.Unbox a => Int -> [(Int, a)] -> SVector a
+fromListSv n ll = SV n $ VU.fromList ll
+
+-- | Dense vectors with unboxed components
+newtype DVector a = DV { dvVec :: VU.Vector a } deriving (Eq, Ord, Generic)
+instance (VU.Unbox a, Serialise a) => Serialise (DVector a)
+instance (VU.Unbox a, Show a) => Show (DVector a) where
+  show (DV vv) = unwords ["DV", show (VU.toList vv)]
+
+fromListDv :: VU.Unbox a => [a] -> DVector a
+fromListDv ll = DV $ VU.fromList ll
+toListDv :: (VU.Unbox a) => DVector a -> [a]
+toListDv (DV v) = VU.toList v
+
+-- | Internal
+--
+-- one projection vector per node (like @annoy@)
+data RT v d a =
+  RBin !d !(v d) !(RT v d a) !(RT v d a)
+  | RTip { _rData :: !a } deriving (Eq, Show, Generic, Functor, Foldable, Traversable)
+makeLensesFor [("_rData", "rData")] ''RT
+instance (NFData (v d), NFData d, NFData a) => NFData (RT v d a)
+
+
+
+-- | Internal
+--
+-- one projection vector per tree level (as suggested in https://www.cs.helsinki.fi/u/ttonteri/pub/bigdata2016.pdf )
+data RPT d a =
+  Bin {
+  _rpThreshold :: !d
+  , _rpMargin :: !(Margin d)
+  , _rpL :: !(RPT d a)
+  , _rpR :: !(RPT d a) }
+  | Tip { _rpData :: a }
+  deriving (Eq, Show, Generic, Functor, Foldable, Traversable)
+instance (Serialise a, Serialise d) => Serialise (RPT d a)
+makeLensesFor [("_rpData", "rpData")] ''RPT
+instance (NFData v, NFData a) => NFData (RPT v a)
+
+-- | Random projection trees
+--
+-- The first type parameter corresponds to a floating point scalar value, the second is the type of the data collected at the leaves of the tree (e.g. lists of vectors)
+--
+-- We keep them separate to leverage the Functor instance for postprocessing and visualization
+--
+-- One projection vector per tree level (as suggested in https://www.cs.helsinki.fi/u/ttonteri/pub/bigdata2016.pdf )
+data RPTree d a = RPTree {
+  _rpVectors :: V.Vector (SVector d) -- ^ one random projection vector per tree level
+  , _rpTree :: RPT d a
+                         } deriving (Eq, Show, Functor, Foldable, Traversable, Generic)
+instance (Serialise d, Serialise a, VU.Unbox d) => Serialise (RPTree d a)
+makeLensesFor [("_rpTree", "rpTree")] ''RPTree
+instance (NFData a, NFData d) => NFData (RPTree d a)
+
+type RPForest d a = IM.IntMap (RPTree d a)
+
+rpTreeData :: Traversal' (RPTree d a) a
+rpTreeData = rpTree . rpData
+
+leaves :: RPTree d a -> [a]
+leaves = (^.. rpTreeData)
+
+-- | Number of tree levels
+levels :: RPTree d a -> Int
+levels (RPTree v _) = VG.length v
+
+-- | Set of data points used to construct the index
+points :: Monoid m => RPTree d m -> m
+points (RPTree _ t) = fold t
+
+-- -- points in 2d
+-- data P a = P !a !a deriving (Eq, Show)
+
+class Scale v where
+  (.*) :: (VU.Unbox a, Num a) => a -> v a -> v a
+instance Scale SVector where
+  a .* (SV n vv) = SV n $ scaleS a vv
+instance Scale VU.Vector where
+  a .* v1 = scaleD a v1
+instance Scale DVector where
+  a .* (DV v1) = DV $ scaleD a v1
+
+-- | Inner product spaces
+--
+-- This typeclass is provided as a convenience for library users to interface their own vector types.
+class (Scale u, Scale v) => Inner u v where
+  inner :: (VU.Unbox a, Num a) => u a -> v a -> a
+  metricL2 :: (VU.Unbox a, Floating a) => u a -> v a -> a
+  (^+^) :: (VU.Unbox a, Num a) => u a -> v a -> u a
+  (^-^) :: (VU.Unbox a, Num a) => u a -> v a -> u a
+
+instance Inner SVector SVector where
+  inner (SV _ v1) (SV _ v2) = innerSS v1 v2
+  metricL2 (SV _ v1) (SV _ v2) = metricSSL2 v1 v2
+  (SV n v1) ^+^ (SV _ v2) = SV n $ sumSS v1 v2
+  (SV n v1) ^-^ (SV _ v2) = SV n $ diffSS v1 v2
+instance Inner SVector VU.Vector where
+  inner (SV _ v1) v2 = innerSD v1 v2
+  metricL2 (SV _ v1) v2 = metricSDL2 v1 v2
+  (SV n v1) ^+^ v2 = SV n $ sumSD v1 v2
+  (SV n v1) ^-^ v2 = SV n $ diffSD v1 v2
+instance Inner SVector DVector where
+  inner (SV _ v1) (DV v2) = innerSD v1 v2
+  metricL2 (SV _ v1) (DV v2) = metricSDL2 v1 v2
+  (SV n v1) ^+^ (DV v2) = SV n $ sumSD v1 v2
+  (SV n v1) ^-^ (DV v2) = SV n $ diffSD v1 v2
+instance Inner DVector DVector where
+  inner (DV v1) (DV v2) = innerDD v1 v2
+  metricL2 (DV v1) (DV v2) = metricDDL2 v1 v2
+  DV v1 ^+^ DV v2 = DV $ VG.zipWith (+) v1 v2
+  DV v1 ^-^ DV v2 = DV $ VG.zipWith (-) v1 v2
+
+(/.) :: (Scale v, VU.Unbox a, Fractional a) => v a -> a -> v a
+v /. a = (1 / a) .* v
+
+normalize :: (VU.Unbox a, Inner v v, Floating a) => v a -> v a
+normalize v = v /. metricL2 v v
+
+
+-- | sparse-sparse inner product
+innerSS :: (VG.Vector u (Int, a), VG.Vector v (Int, a), VU.Unbox a, Num a) =>
+           u (Int, a) -> v (Int, a) -> a
+innerSS vv1 vv2 = go 0 0
+  where
+    nz1 = VG.length vv1
+    nz2 = VG.length vv2
+    go i1 i2
+      | i1 >= nz1 || i2 >= nz2 = 0
+      | otherwise =
+          let
+            (il, xl) = vv1 VG.! i1
+            (ir, xr) = vv2 VG.! i2
+          in case il `compare` ir of
+            EQ -> (xl * xr +) $ go (succ i1) (succ i2)
+            LT -> go (succ i1) i2
+            GT -> go i1 (succ i2)
+
+-- | sparse-dense inner product
+innerSD :: (Num a, VG.Vector u (Int, a), VG.Vector v a, VU.Unbox a) =>
+           u (Int, a) -> v a -> a
+innerSD vv1 vv2 = go 0
+  where
+    nz1 = VG.length vv1
+    nz2 = VG.length vv2
+    go i
+      | i >= nz1 || i >= nz2 = 0
+      | otherwise =
+          let
+            (il, xl) = vv1 VG.! i
+            xr       = vv2 VG.! il
+          in
+            (xl * xr +) $ go (succ i)
+
+innerDD :: (VG.Vector v a, Num a) => v a -> v a -> a
+innerDD v1 v2 = VG.sum $ VG.zipWith (*) v1 v2
+
+
+-- | Vector distance induced by the L2 norm (sparse-sparse)
+metricSSL2 :: (Floating a, VG.Vector u a, VU.Unbox a, VG.Vector u (Int, a), VG.Vector v (Int, a)) =>
+              u (Int, a) -> v (Int, a) -> a
+metricSSL2 u v = sqrt $ VG.sum $ VG.map (\(_, x) -> x ** 2) duv
+  where
+    duv = u `diffSS` v
+
+-- | Vector distance induced by the L2 norm (sparse-dense)
+metricSDL2 :: (Floating a, VG.Vector v1 a, VU.Unbox a,
+                VG.Vector v1 (Int, a), VG.Vector v2 a) =>
+              v1 (Int, a) -> v2 a -> a
+metricSDL2 u v = sqrt $ VG.sum $ VG.map (\(_, x) -> x ** 2) duv
+  where
+    duv = u `diffSD` v
+
+-- | Vector distance induced by the L2 norm (dense-dense)
+metricDDL2 :: (Floating a, VG.Vector v a) => v a -> v a -> a
+metricDDL2 u v = sqrt $ VG.sum $ VG.map (** 2) duv
+  where
+    duv = VG.zipWith (-) u v
+
+scaleD :: (VG.Vector v b, Num b) => b -> v b -> v b
+scaleD a vv = VG.map (* a) vv
+
+scaleS :: (VG.Vector v (a, b), Num b) => b -> v (a, b) -> v (a, b)
+scaleS a vv = VG.map (\(i, x) -> (i, a * x)) vv
+
+-- | Vector sum
+sumSD :: (VG.Vector u (Int, a), VG.Vector v a, VU.Unbox a, Num a) =>
+         u (Int, a) -> v a -> u (Int, a)
+sumSD = binSD (-)
+
+-- | Vector sum
+sumSS :: (VG.Vector u (Int, a), VG.Vector v (Int, a), VU.Unbox a, Num a) =>
+         u (Int, a) -> v (Int, a) -> u (Int, a)
+sumSS = binSS (+) 0 
+
+-- | Vector difference
+diffSD :: (VG.Vector u (Int, a), VG.Vector v a, VU.Unbox a, Num a) =>
+          u (Int, a) -> v a -> u (Int, a)
+diffSD = binSD (-)
+
+-- | Vector difference
+diffSS :: (VG.Vector u (Int, a), VG.Vector v (Int, a), VU.Unbox a, Num a) =>
+          u (Int, a) -> v (Int, a) -> u (Int, a)
+diffSS = binSS (-) 0
+
+-- | Binary operation on 'SVector' s
+binSS :: (VG.Vector u (Int, a), VG.Vector v (Int, a), VU.Unbox a) =>
+         (a -> a -> a) -> a -> u (Int, a) -> v (Int, a) -> u (Int, a)
+binSS f z vv1 vv2 = VG.unfoldr go (0, 0)
+  where
+    nz1 = VG.length vv1
+    nz2 = VG.length vv2
+    go (i1, i2)
+      | i1 >= nz1 || i2 >= nz2 = Nothing
+      | otherwise =
+          let
+            (il, xl) = vv1 VG.! i1
+            (ir, xr) = vv2 VG.! i2
+          in case il `compare` ir of
+            EQ -> Just ((il, f xl xr), (succ i1, succ i2))
+            LT -> Just ((il, f xl z ), (succ i1, i2     ))
+            GT -> Just ((ir, f z  xr), (i1     , succ i2))
+
+
+
+binSD :: (VG.Vector u (Int, a), VG.Vector v a, VU.Unbox a) =>
+         (a -> a -> a) -> u (Int, a) -> v a -> u (Int, a)
+binSD f vv1 vv2 = VG.unfoldr go 0
+  where
+    nz1 = VG.length vv1
+    nz2 = VG.length vv2
+    go i
+      | i >= nz1 || i >= nz2 = Nothing
+      | otherwise = Just ((il, y), succ i)
+          where
+            (il, xl) = vv1 VG.! i
+            xr       = vv2 VG.! il
+            y = f xl xr
+
+
+-- | Partition the data wrt the median value of the inner product
+partitionAtMedian :: (Ord a, Inner u v, VU.Unbox a, Fractional a) =>
+                     u a -- ^ projection vector
+                  -> V.Vector (v a) -- ^ dataset (3 or more elements)
+                  -> (a, Margin a, V.Vector (v a), V.Vector (v a)) -- ^ median, margin, smaller, larger
+partitionAtMedian r xs = (thr, margin, ll, rr)
+  where
+    (ll, rr) = (VG.take nh xs', VG.drop nh xs')
+    -- (pjl, pjr) = (VG.head inns, VG.last inns) -- (min, max) inner product values
+    (mgl, mgr) = (inns VG.! (nh - 1), inns VG.! (nh + 1))
+    margin = Margin (Max mgl) (Min mgr)
+    -- marginL = mgl / (pjr - pjl) -- lower bound of margin, normalized to range
+    -- marginR = mgr / (pjr - pjl) -- upper bound of margin, normalized to range
+    thr = inns VG.! nh -- inner product threshold
+    n = VG.length xs -- total data size
+    nh = n `div` 2 -- size of left partition
+    projs = sortByVG snd $ VG.map (\x -> (x, r `inner` x)) xs
+    (xs', inns) = VG.unzip projs
+
+sortByVG :: (VG.Vector v a, Ord b) => (a -> b) -> v a -> v a
+sortByVG f v = runST $ do
+  vm <- VG.thaw v
+  V.sortBy (comparing f) vm
+  VG.freeze vm
+
+
+
+
+
+
+
+-- data Avg a = Avg {
+--   avgCount :: !(Sum Int)
+--   , avgSum :: !(Sum a)
+--                  }
+-- average :: (Foldable t, Fractional a) => t a -> a
+-- average = getAvg . foldl' bumpAvg mempty
+-- bumpAvg :: Num a => Avg a -> a -> Avg a
+-- bumpAvg aa x = Avg (Sum 1) (Sum x) <> aa
+-- instance (Num a) => Semigroup (Avg a) where
+--   Avg c0 s0 <> Avg c1 s1 = Avg (c0<>c1) (s0<>s1)
+-- instance (Num a) => Monoid (Avg a) where
+--   mempty = Avg mempty mempty
+-- getAvg :: Fractional a => Avg a -> a
+-- getAvg (Avg c s) = getSum s / fromIntegral (getSum c)
+
+
+-- -- | Label a value with a unique identifier
+-- -- labelId
+-- newtype LabelT m a = LabelT {unLabelT :: StateT Integer m a} deriving (Functor, Applicative, Monad, MonadState Integer, MonadIO)
+-- type Label = LabelT Identity
+-- runLabelT :: (Monad m) => LabelT m a -> m a
+-- runLabelT = flip evalStateT 0 . unLabelT
+-- label :: Monad m => a -> LabelT m (Id a)
+-- label x = LabelT $ do { i <- get ; put (i + 1); pure (Id x i)}
+-- data Id a = Id { _idD :: a , _idL :: !Integer } deriving (Eq, Show, Functor, Foldable, Traversable, Generic)
+-- instance NFData a => NFData (Id a)
+-- makeLensesFor [("_idD", "idD")] ''Id
+-- instance (Eq a) => Ord (Id a) where
+--   Id _ u1 <= Id _ u2 = u1 <= u2
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,1 @@
+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}