rp-tree 0.1.0.0 → 0.2.0.0

raw patch · 10 files changed

+654/−142 lines, 10 filesdep +benchpressdep +mnist-idx-conduitdep +splitmixdep ~basedep ~splitmix-distributionsdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: benchpress, mnist-idx-conduit, splitmix

Dependency ranges changed: base, splitmix-distributions, vector

API changes (from Hackage documentation)

Files

• Changelog.md +6/−0
• README.md +2/−2
• app/Main.hs +32/−13
• bench/time/Main.hs +277/−0
• rp-tree.cabal +29/−4
• src/Data/RPTree.hs +52/−39
• src/Data/RPTree/Conduit.hs +107/−68
• src/Data/RPTree/Gen.hs +42/−1
• src/Data/RPTree/Internal.hs +57/−15
• src/Data/RPTree/Internal/Testing.hs +50/−0

@@ -0,0 +1,6 @@+0.2++- fix 'candidates' such that 'knn' now does the right thing+- now 'knn' accepts a distance function as parameter as well+- add I/O functionality+- some time benchmarks

@@ -1,5 +1,5 @@ # rp-tree -![scatterplot](/r/scatter.png "scatterplot")+![rp-tree](https://github.com/ocramz/rp-tree/blob/main/r/scatter.png ) -Random projection trees+Random projection trees for approximate nearest neighbor search in high-dimensional vector spaces

@@ -3,9 +3,12 @@ {-# options_ghc -Wno-unused-imports #-} module Main where +import Control.Monad (replicateM) import Data.Foldable (fold)+import Data.Functor (void)+ -- conduit-import qualified Data.Conduit as C (ConduitT, runConduit, yield, await)+import qualified Data.Conduit as C (ConduitT, runConduit, yield, await, transPipe) 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)@@ -13,25 +16,46 @@ import qualified Data.IntMap as IM (IntMap, fromList, insert, lookup, map, mapWithKey, traverseWithKey, foldlWithKey, foldrWithKey) -- exceptions import Control.Monad.Catch (MonadThrow(..))+-- mnist-idx-conduit+import Data.IDX.Conduit (sourceIdxSparse, sBufSize, sNzComponents) -- splitmix-distributions import System.Random.SplitMix.Distributions (Gen, GenT, sample, sampleT, bernoulli, normal) -- transformers import Control.Monad.Trans.State.Lazy (State, get, put, evalState)+import Control.Monad.Trans.Class (MonadTrans(..)) -- 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)+import Data.RPTree (knn, candidates, Inner(..), RPTree, RPForest, leaves, SVector, fromListSv, DVector, fromListDv, dense, writeCsv, tree, forest, dataSource, sparse, normal2, normalSparse2)+import Data.RPTree.Internal.Testing (datS, datD)  main :: IO ()-main = do+main = do -- putStrLn "hello!"   let-    n = 10000-  -- renderTree0 n-  -- renderTree1 n-  undefined -- FIXME+    n = 1000+    maxd = 3+    minl = 10+    ntree = 10+    d = 100+    pnz = 0.3+    chunk = 20+    src = datS n d pnz+    seed = 1234+  (q, tts) <- sampleT seed $ do+    tts <- C.runConduit $+           forest seed maxd minl ntree chunk pnz d (liftC src)+    q <- sparse 0.3 d (normal 0.1 0.6)+    pure (q, tts)+  let+    res = knn metricL2 1 tts q+  print res  +liftC = C.transPipe lift+++ -- renderTree0 :: Int -> IO () renderTree0 tt = do   let csvrows = V.toList $ fold $ flip evalState A $ traverse labeledV tt -- (tree0 n)@@ -115,9 +139,4 @@ 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+

@@ -0,0 +1,277 @@+{-# LANGUAGE FlexibleContexts #-}+{-# options_ghc -Wno-unused-imports #-}+module Main where++import Control.Exception (bracket)+import Control.Monad (forM, forM_)+import Control.Monad.IO.Class (MonadIO(..))+import GHC.Word (Word8, Word64)+import System.CPUTime (getCPUTime)++-- -- benchpress+-- import Test.BenchPress (Stats(..), benchmark, printDetailedStats)+-- conduit+import Conduit (runResourceT, MonadResource)+import qualified Data.Conduit as C (ConduitT, runConduit, runConduitRes, yield, await, transPipe)+import Data.Conduit ((.|))+import qualified Data.Conduit.Combinators as C (map, mapM, scanl, scanlM, last, print, takeExactly)+-- deepseq+import Control.DeepSeq (NFData(..), force)+-- exceptions+import Control.Monad.Catch (MonadThrow(..))+-- mnist-idx-conduit+import Data.IDX.Conduit (sourceIdxSparse, sBufSize, sNzComponents)+-- splitmix-distributions+import System.Random.SplitMix.Distributions (GenT, sampleT, sample, samples)++-- mtl+import Control.Monad.Trans.Class (MonadTrans(..))++-- vector+import qualified Data.Vector as V (Vector, replicateM, fromList)+import qualified Data.Vector.Unboxed as VU (Unbox, Vector, map)++import Data.RPTree (tree, forest, recallWith, knn, fromVectorSv, fromListSv, RPForest, RPTree, SVector, Inner(..), normalSparse2, liftC)+import Data.RPTree.Internal.Testing (BenchConfig(..), randSeed, datD, datS)++main :: IO ()+main = do -- putStrLn "hello!"+  binMixFQBench+  -- mnistBench++benchConfigs :: String -- ^ description of the experiment+             -> Int -- ^ dimension of the projection vectors+             -> [BenchConfig]+benchConfigs descr pdim = [ BenchConfig descr maxd minl nt chunk nzd pdim n nq+                          | maxd <- [5],+                            minl <- [10],+                            nt <- [3],+                            chunk <- [100],+                            nzd <- [0.2],+                            n <- [1000],+                            nq <- [10, 100]+                          ]++-- -- Binary mixture++binMixFQBench :: IO ()+binMixFQBench = do+  let+    cfgs = benchConfigs "binary mixture of sparse Gaussian RVs" 1000+  forM_ cfgs $ \cfg -> do+    stats <- binMixFQBench1 cfg+    print cfg+    -- printDetailedStats stats+    print stats+    pure stats+  -- print s++-- | Measure recall @ 10 and mean time+binMixFQBench1 :: BenchConfig -> IO (Double, Double)+binMixFQBench1 cfg = forestBenchGen seed (datS n d nzData) act 2 cfg+  where+    n = bcDataSize cfg+    d = bcVectorDim cfg+    nq = bcNumQueryPoints cfg+    -- pnz = bcNZDensity cfg -- nz density of proj vectors+    nzData = 0.8 -- nz density of data +    k = 10 -- number of ANN's to return+    seed = 1234+    qs = samples nq seed $ normalSparse2 nzData d+    act tt = do+      -- pure $! recallWith metricL2 tt k `map` qs+      let+        recs = recallWith metricL2 tt k `map` qs+        r = mean recs+      pure $! r++++-- -- MNIST dataset++mnistBench :: IO ()+mnistBench = do+  let+    cfgs = benchConfigs "MNIST dataset" 784+    mnfpath = "assets/mnist/train-images-idx3-ubyte"+  forM_ cfgs $ \cfg -> do+    stats <- mnistFQBench1 mnfpath cfg+    print cfg+    print stats++-- | Measure recall @ 10 and mean time+mnistFQBench1 :: FilePath -> BenchConfig -> IO (Double, Double)+mnistFQBench1 fp cfg = forestBench (mnist fp n) act 1 cfg+  where+    n = bcDataSize cfg+    k = 10 -- number of ANN's to return+    d = bcVectorDim cfg+    nzData = 0.8 -- nz density of data +    act x = do+      tt <- runResourceT x+      let q = sample 1234 $ normalSparse2 nzData d+      pure $! recallWith metricL2 tt k q++mnist :: MonadResource m =>+         FilePath -- path to uncompressed MNIST IDX data file+      -> Int -- number of data items+      -> C.ConduitT a (SVector Double) m ()+mnist fp n = C.takeExactly n src+  where+    src = sourceIdxSparse fp .|+          C.map (\r -> fromVectorSv (sBufSize r) (VU.map f $ sNzComponents r))+    f (i, x) = (i, toUnitRange x)++toUnitRange :: Word8 -> Double+toUnitRange w8 = fromIntegral w8 / 255+++++-- -- UTILS++-- | runs a benchmark on a newly created RPForest initialized with a random seed+forestBench :: (MonadThrow m, Inner SVector v) =>+               C.ConduitT () (v Double) m ()+            -> (m (RPForest Double (V.Vector (v Double))) -> IO c) -- ^ allows for both deterministic and random data sources+            -> Int  -- ^ number of replicates+            -> BenchConfig+            -> IO (c, Double) -- ^ result, mean wall-clock time measurement+forestBench src go n cfg = benchmark n setup (const $ pure ()) go+  where+    setup = do+      s <- randSeed+      -- let src' = C.transPipe lift src+      pure $ growForest s cfg src+++-- -- forestBenchRes :: (MonadResource m, Inner SVector v) =>+-- --                   C.ConduitT+-- --                   ()+-- --                   (v Double)+-- --                   m+-- --                   ()+-- --                -> (RPForest Double (V.Vector (v Double)) -> IO c)+-- --                -> Int+-- --                -> BenchConfig+-- --                -> IO (c, Double)+forestBenchRes src go n cfg = benchmarkM n setup go+  where+    setup = do+      s <- randSeed+      runResourceT $ growForest s cfg src++forestBenchGen :: (MonadIO m, Inner SVector v, NFData (v Double)) =>+                  Word64+               -> C.ConduitT () (v Double) (GenT m) ()+               -> (RPForest Double (V.Vector (v Double)) -> m a2)+               -> Int+               -> BenchConfig+               -> m (a2, Double)+forestBenchGen seed src go n cfg = benchmarkM n setup go+  where+    setup = do+      s <- randSeed+      x <- sampleT seed $ growForest s cfg src+      pure $ force x++++treeBench :: (Monad m, Inner SVector v) =>+             C.ConduitT () (v Double) m ()+          -> (m (RPTree Double (V.Vector (v Double))) -> IO c)+          -> Int+          -> BenchConfig+          -> IO (c, Double)+treeBench src go n cfg = benchmark n setup (const $ pure ()) go+      where+        setup = do+          s <- randSeed+          -- let src' = C.transPipe lift src+          pure $ growTree s cfg src++growTree :: (Monad m, Inner SVector v) =>+            Word64+         -> BenchConfig+         -> C.ConduitT () (v Double) m ()+         -> m (RPTree Double (V.Vector (v Double)))+growTree seed (BenchConfig _ maxd minl _ chunksize pnz pdim _ _) =+  tree seed maxd minl chunksize pnz pdim++growForest :: (Monad m, Inner SVector v) =>+              Word64+           -> BenchConfig+           -> C.ConduitT () (v Double) m ()+           -> m (RPForest Double (V.Vector (v Double)))+growForest seed (BenchConfig _ maxd minl ntrees chunksize pnz pdim _ _) =+  forest seed maxd minl ntrees chunksize pnz pdim++-- growForest' seed (BenchConfig _ maxd minl ntrees chunksize pnz pdim _) =+--   forest' seed maxd minl ntrees chunksize pnz pdim+++-- -- adapted from 'benchpress', until https://github.com/WillSewell/benchpress/issues/9 is merged+benchmark :: Int -> IO a -> (a -> IO b) -> (a -> IO c) -> IO (c, Double)+benchmark iters setup teardown action =+  if iters < 1+    then error "benchmark: iters must be greater than 0"+    else do+      (vals, cpuTimes) <- unzip `fmap` go iters+      let tcm = mean cpuTimes+          v = head vals+      return (v, tcm)+      where+        go 0 = return []+        go n = do+          elapsed <- bracket setup teardown $ \a -> do+            startCpu <- getCPUTime+            x <- action a+            endCpu <- getCPUTime+            return (x+                   ,picosToMillis $! endCpu - startCpu)+          timings <- go $! n - 1+          return $ elapsed : timings++benchmarkM :: (MonadIO m) =>+              Int -> m t -> (t -> m a2) -> m (a2, Double)+benchmarkM iters setup action =+  if iters < 1+    then error "benchmark: iters must be greater than 0"+    else do+      (vals, cpuTimes) <- unzip `fmap` go iters+      let tcm = mean cpuTimes+          v = head vals+      return (v, tcm)+      where+        go 0 = return []+        go n = do+          a <- setup+          elapsed <- do+            startCpu <- liftIO getCPUTime+            x <- action a+            endCpu <- liftIO getCPUTime+            return (x+                   ,picosToMillis $! endCpu - startCpu)+          timings <- go $! n - 1+          return $ elapsed : timings+++-- | Converts picoseconds to milliseconds.+picosToMillis :: Integer -> Double+picosToMillis t = realToFrac t / (10^(9 :: Int))++-- | Numerically stable mean.+mean :: Floating a => [a] -> a+mean = go 0 0+    where+      go :: Floating a => a -> Int -> [a] -> a+      go m _ []     = m+      go m n (x:xs) = go (m + (x - m) / fromIntegral (n + 1)) (n + 1) xs+++++mnistV0 :: SVector Double+mnistV0 = fromListSv 784 (map (\(i, x) -> (i, toUnitRange x)) cs)+  where+    cs = [(208,55),(209,148),(210,210),(211,253),(212,253),(213,113),(214,87),(215,148),(216,55),(235,87),(236,232),(237,252),(238,253),(239,189),(240,210),(241,252),(242,252),(243,253),(244,168),(261,4),(262,57),(263,242),(264,252),(265,190),(266,65),(267,5),(268,12),(269,182),(270,252),(271,253),(272,116),(289,96),(290,252),(291,252),(292,183),(293,14),(296,92),(297,252),(298,252),(299,225),(300,21),(316,132),(317,253),(318,252),(319,146),(320,14),(324,215),(325,252),(326,252),(327,79),(343,126),(344,253),(345,247),(346,176),(347,9),(350,8),(351,78),(352,245),(353,253),(354,129),(370,16),(371,232),(372,252),(373,176),(377,36),(378,201),(379,252),(380,252),(381,169),(382,11),(398,22),(399,252),(400,252),(401,30),(402,22),(403,119),(404,197),(405,241),(406,253),(407,252),(408,251),(409,77),(426,16),(427,231),(428,252),(429,253),(430,252),(431,252),(432,252),(433,226),(434,227),(435,252),(436,231),(455,55),(456,235),(457,253),(458,217),(459,138),(460,42),(461,24),(462,192),(463,252),(464,143),(489,62),(490,255),(491,253),(492,109),(517,71),(518,253),(519,252),(520,21),(546,253),(547,252),(548,21),(573,71),(574,253),(575,252),(576,21),(601,106),(602,253),(603,252),(604,21),(629,45),(630,255),(631,253),(632,21),(658,218),(659,252),(660,56),(686,96),(687,252),(688,189),(689,42),(714,14),(715,184),(716,252),(717,170),(718,11),(743,14),(744,147),(745,252),(746,42)]

@@ -1,5 +1,5 @@ name:                rp-tree-version:             0.1.0.0+version:             0.2.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@@ -11,6 +11,7 @@ category:            Data Mining, Data Structures, Machine Learning, Data build-type:          Simple extra-source-files:  README.md+                     Changelog.md cabal-version:       >=1.10 tested-with:         GHC == 8.10.4 @@ -19,7 +20,7 @@   ghc-options:         -Wall   hs-source-dirs:      src   exposed-modules:     Data.RPTree-+                       Data.RPTree.Internal.Testing   other-modules:       Data.RPTree.Internal                        Data.RPTree.Gen                        Data.RPTree.Draw@@ -41,8 +42,11 @@                      -- , ulid                      , vector >= 0.12.1.2                      , vector-algorithms-                       -- debug+                     -- -- DEBUG+                     , benchpress                      , hspec+                     , mnist-idx-conduit+                     , splitmix >= 0.1.0.3  test-suite spec   default-language:    Haskell2010@@ -56,16 +60,37 @@                      , QuickCheck                      , splitmix-distributions +benchmark bench-time+  default-language:    Haskell2010+  ghc-options:         -threaded -O2+  type:                exitcode-stdio-1.0+  hs-source-dirs:      bench/time+  main-is:             Main.hs+  build-depends:       base+                     , benchpress+                     , conduit+                     , deepseq >= 1.4.4.0+                     , exceptions+                     , mnist-idx-conduit+                     , mtl+                     , rp-tree+                     , splitmix >= 0.1.0.3+                     , splitmix-distributions+                     , transformers+                     , vector+ executable rp-tree   default-language:    Haskell2010-  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  ghc-options:         -threaded -O2   hs-source-dirs:      app   main-is:             Main.hs   build-depends:       base                      , conduit                      , containers                      , exceptions+                     , mnist-idx-conduit                      , rp-tree+                     , splitmix                      , splitmix-distributions                      , transformers                      , vector

@@ -15,40 +15,45 @@ -} module Data.RPTree (   -- * Construction-  forest+  tree, forest   -- * Query   , knn-  -- , nearest+  -- * I/O+  , serialiseRPForest+  , deserialiseRPForest   -- * Validation-  , recall+  , recallWith   -- * Access   , levels, points, leaves, candidates   -- * Types   -- ** RPTree   , RPTree, RPForest-  -- -- *** internal-  -- , RPT-  -- -- ** RT-  -- , RT   -- *-  , SVector, fromListSv-  , DVector, fromListDv-  -- * inner product+  , SVector, fromListSv, fromVectorSv+  , DVector, fromListDv, fromVectorDv+  -- * Vector space types   , Inner(..), Scale(..)-  --   -- ** helpers for implementing Inner instances-  --   -- *** inner product-  -- , innerSS, innerSD-  --   -- *** L2 distance-  -- , metricSSL2, metricSDL2+    -- ** helpers for implementing Inner instances+    -- *** inner product+  , innerSS, innerSD, innerDD+    -- *** L2 distance+  , metricSSL2, metricSDL2+  -- *** Scale+  , scaleS, scaleD   -- * Conduit   , dataSource   -- * Random generation   -- ** vector   , sparse, dense+  , normal2+   -- * Rendering   , draw   -- * CSV   , writeCsv+  -- * Testing+  , randSeed, BenchConfig(..), normalSparse2+  , liftC   ) where  import Control.Monad (replicateM)@@ -83,10 +88,10 @@ -- 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.Conduit (tree, forest, dataSource, liftC)+import Data.RPTree.Gen (sparse, dense, normal2, normalSparse2)+import Data.RPTree.Internal (RPTree(..), RPForest, RPT(..), levels, points, leaves, RT(..), Inner(..), Scale(..), scaleS, scaleD, (/.), innerDD, innerSD, innerSS, metricSSL2, metricSDL2, SVector(..), fromListSv, fromVectorSv, DVector(..), fromListDv, fromVectorDv, partitionAtMedian, Margin, getMargin, sortByVG, serialiseRPForest, deserialiseRPForest)+import Data.RPTree.Internal.Testing (BenchConfig(..), randSeed) import Data.RPTree.Draw (draw, writeCsv)  @@ -103,32 +108,37 @@   -- | 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+recallWith :: (Inner SVector v, VU.Unbox a, Fractional a, Ord a, Ord d, Ord p) =>+              (p -> v a -> d)+           -> RPForest a (V.Vector p)+           -> Int -- ^ k : number of nearest neighbors to consider+           -> v a -- ^ query point+           -> a+recallWith distf tt k q = sum rs / fromIntegral n   where-    rs = fmap (\t -> recall1 t k q) tt+    rs = fmap (\t -> recallWith1 distf 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+-- -- | Recall computed wrt the Euclidean distance+-- recallEuclid :: (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+-- recallEuclid = recallWith metricL2 -recallWith :: (Fractional a1, Inner SVector v, Ord d, VU.Unbox d,+recallWith1 :: (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+              (a2 -> v d -> a3) -- ^ distance function+           -> RPTree d (V.Vector a2)+           -> Int -- ^ k : number of nearest neighbors to consider+           -> v d -- ^ query point+           -> a1+recallWith1 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+    kk = S.fromList $ map fst $ take k dists -- first k nn's     aa = set $ candidates tt q     aintk = aa `S.intersection` kk @@ -137,6 +147,7 @@   +{-# SCC candidates #-} -- | Retrieve points nearest to the query -- -- in case of a narrow margin, collect both branches of the tree@@ -153,11 +164,13 @@         r = rvs VG.! ixLev         proj = r `inner` x         i' = succ ixLev+        dl = abs (mglo - proj) -- left margin+        dr = abs (mghi - proj) -- right margin       if | proj < thr &&-           mglo > mghi -> go i' ltree <> go i' rtree+           dl > dr -> go i' ltree <> go i' rtree          | proj < thr  -> go i' ltree          | proj > thr &&-           mglo < mghi -> go i' ltree <> go i' rtree+           dl < dr -> go i' ltree <> go i' rtree          | otherwise   -> go i' rtree  

@@ -1,23 +1,28 @@ {-# language DeriveDataTypeable #-} {-# LANGUAGE FlexibleContexts #-}+{-# language BangPatterns #-} {-# options_ghc -Wno-unused-imports #-} {-# options_ghc -Wno-unused-top-binds #-}-module Data.RPTree.Conduit (+module Data.RPTree.Conduit+  (tree,   forest   -- ** helpers   , dataSource-  ) where+  , liftC+  )+where  import Control.Monad (replicateM)+import Data.Functor (void) import GHC.Word (Word64)  -- conduit-import qualified Data.Conduit as C (ConduitT, runConduit, yield, await)+import qualified Data.Conduit as C (ConduitT, runConduit, yield, await, transPipe) 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)+import qualified Data.Conduit.Combinators as C (map, mapM, last, scanl, print, foldl)+import qualified Data.Conduit.List as C (chunksOf, unfold, unfoldM, mapAccum) -- containers-import qualified Data.IntMap as IM (IntMap, fromList, insert, lookup, map, mapWithKey, traverseWithKey, foldlWithKey, foldrWithKey, intersectionWith)+import qualified Data.IntMap.Strict as IM (IntMap, fromList, insert, lookup, map, mapWithKey, traverseWithKey, foldlWithKey, foldrWithKey, intersectionWith) -- exceptions import Control.Monad.Catch (MonadThrow(..)) -- mtl@@ -37,19 +42,8 @@ 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)+liftC :: (Monad m, MonadTrans t) => C.ConduitT i o m r -> C.ConduitT i o (t m) r+liftC = C.transPipe lift  -- | Populate a tree from a data stream --@@ -62,33 +56,34 @@ -- * 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 :: (Monad 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"+  t <- C.runConduit $ src .|+                      insertC maxDepth minLeaf n rvs+  pure $ RPTree rvs t ++++ -- | 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)+        -> C.ConduitT (v d) o m (RPT d (V.Vector (v d))) +insertC maxDepth minLeaf n rvs = chunkedAccum n z (insert maxDepth minLeaf rvs)   where     z = Tip mempty @@ -105,35 +100,30 @@ -- * 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 :: (Monad 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"+  ts <- C.runConduit $ src .|+                       insertMultiC maxd minl chunksize rvss+  pure $ IM.intersectionWith RPTree rvss ts     + 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@@ -141,15 +131,16 @@              -> IM.IntMap (v1 (u d)) -- one entry per tree              -> C.ConduitT                 (v d)-                (IM.IntMap (RPT d (V.Vector (v d))))+                o                 m-                ()-insertMultiC maxd minl n rvss = chunked n im0 (insertMulti maxd minl rvss)+                (IM.IntMap (RPT d (V.Vector (v d))))+insertMultiC maxd minl n rvss = chunkedAccum n im0 (insertMulti maxd minl rvss)   where     im0 = IM.map (const z) rvss     z = Tip mempty  +{-# SCC insertMulti #-} insertMulti :: (Ord d, Inner u v, VU.Unbox d, Fractional d, VG.Vector v1 (u d)) =>                Int             -> Int@@ -158,10 +149,11 @@             -> 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+  flip IM.mapWithKey tacc $ \ !i !t -> case IM.lookup i rvss of+                                      Just !rvs -> insert maxd minl rvs t xs                                       _        -> t +{-# SCC insert #-} insert :: (VG.Vector v1 (u d), Ord d, Inner u v, VU.Unbox d, Fractional d) =>           Int -- ^ max tree depth        -> Int -- ^ min leaf size@@ -172,7 +164,7 @@ insert maxDepth minLeaf rvs = loop 0   where     z = Tip mempty-    loop ixLev tt xs =+    loop ixLev !tt xs =       let         r = rvs VG.! ixLev       in@@ -203,13 +195,60 @@                 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 .|+-- | Aggregate the input stream in chunks of a given size (semantics of 'C.chunksOf'), and fold over the resulting stream building up an accumulator structure (e.g. a tree)+chunkedAccum :: (Monad m) =>+                Int -- ^ chunk size+             -> t -- ^ initial accumulator state+             -> (t -> V.Vector a -> t)+             -> C.ConduitT a o m t+chunkedAccum n z f = C.chunksOf n .|                      C.map V.fromList .|-                     C.scanl f z -- .|+                     C.foldl f z +-- | 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) ++++-- -- sinks++-- tree' :: (Monad 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) o m (RPTree Double (V.Vector (v Double)))+-- tree' seed maxDepth minLeaf n pnz dim = do+--   let+--     rvs = sample seed $ V.replicateM maxDepth (sparse pnz dim stdNormal)+--   t <- insertC maxDepth minLeaf n rvs+--   pure $ RPTree rvs t++-- forest' :: (Monad 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) o m (IM.IntMap (RPTree Double (V.Vector (v Double))))+-- forest' seed maxd minl ntrees chunksize pnz dim = do+--   let+--     rvss = sample seed $ do+--       rvs <- replicateM ntrees $ V.replicateM maxd (sparse pnz dim stdNormal)+--       pure $ IM.fromList $ zip [0 .. ] rvs+--   ts <- insertMultiC maxd minl chunksize rvss+--   pure $ IM.intersectionWith RPTree rvss ts

@@ -13,7 +13,7 @@ import Control.Monad.Trans.Class (MonadTrans(..)) import Control.Monad.State (MonadState(..), modify) -- splitmix-distribitions-import System.Random.SplitMix.Distributions (Gen, GenT, stdUniform, bernoulli)+import System.Random.SplitMix.Distributions (Gen, GenT, stdUniform, bernoulli, exponential, normal, discrete, categorical) -- transformers import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState) -- vector@@ -96,6 +96,47 @@   let ix = floor (fromIntegral k * u)   pure $ IM.insert ix y imm +++++++-- mixtures++mixtureN :: Monad m => [(Double, GenT m b)] -> GenT m b+mixtureN pgs = go+  where+    (ps, gs) = unzip pgs+    go = do+      miix <- categorical ps+      case miix of+        Nothing -> gs !! 0+        Just i  -> do+          let p = gs !! i+          p+++normalSparse2 :: Monad m => Double -> Int -> GenT m (SVector Double)+normalSparse2 pnz d = do+  b <- bernoulli 0.5+  if b+    then sparse pnz d (normal 0 0.5)+    else sparse pnz d (normal 2 0.5)++normalDense2 :: Monad m => Int -> GenT m (DVector Double)+normalDense2 d = do+  b <- bernoulli 0.5+  if b+    then dense d (normal 0 0.5)+    else dense d (normal 2 0.5)++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   -- | Generate a sparse random vector with a given nonzero density and components sampled from the supplied random generator

@@ -10,6 +10,7 @@ {-# language MultiParamTypeClasses #-} {-# language GeneralizedNewtypeDeriving #-} {-# language TemplateHaskell #-}+{-# LANGUAGE BangPatterns        #-} {-# options_ghc -Wno-unused-imports #-} module Data.RPTree.Internal where @@ -17,7 +18,7 @@ import Control.Monad.IO.Class (MonadIO(..)) import Control.Monad.ST (runST) import Data.Function ((&))-import Data.Foldable (fold, foldl')+import Data.Foldable (fold, foldl', toList) import Data.Functor.Identity (Identity(..)) import Data.List (nub) import Data.Monoid (Sum(..))@@ -26,8 +27,10 @@ import Data.Typeable (Typeable) import GHC.Generics (Generic) +-- bytestring+import qualified Data.ByteString.Lazy as LBS (ByteString, toStrict, fromStrict) -- containers-import qualified Data.IntMap as IM (IntMap)+import qualified Data.IntMap.Strict as IM (IntMap, fromList) -- deepseq import Control.DeepSeq (NFData(..)) -- microlens@@ -36,7 +39,7 @@ -- mtl import Control.Monad.State (MonadState(..), modify) -- serialise-import Codec.Serialise (Serialise(..))+import Codec.Serialise (Serialise(..), serialise, deserialiseOrFail) -- transformers import Control.Monad.Trans.State (StateT(..), runStateT, evalStateT, State, runState, evalState) -- vector@@ -47,6 +50,13 @@ -- vector-algorithms import qualified Data.Vector.Algorithms.Merge as V (sortBy) +-- | Pair a datum with a vector embedding+data Embed v e a = Embed {+  eEmbed :: !(v e) -- ^ the embedding is a vector+  , eData :: !a+                       } deriving (Eq, Ord, Show, Generic, Functor)+instance (NFData (v e), NFData a) => NFData (Embed v e a)+instance (Serialise (v e), Serialise a) => Serialise (Embed v e a)  -- | Exceptions data RPTError =@@ -59,8 +69,8 @@  -- | Bounds around the cutting plane data Margin a = Margin {-  cMarginLow :: Max a -- ^ lower bound on the cut point-  , cMarginHigh :: Min a -- ^ upper bound+  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)@@ -72,7 +82,8 @@   -- | Sparse vectors with unboxed components-data SVector a = SV { svDim :: !Int, svVec :: VU.Vector (Int, a) } deriving (Eq, Ord, Generic)+data SVector a = SV { svDim :: {-# UNPACK #-} !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)]@@ -80,15 +91,29 @@  fromListSv :: VU.Unbox a => Int -> [(Int, a)] -> SVector a fromListSv n ll = SV n $ VU.fromList ll+-- | (Unsafe) Pack a 'SVector' from its vector dimension and components+--+-- Note : the relevant invariants are not checked :+--+-- * vector components are _assumed_ to be in increasing order+--+-- * vector dimension is larger than any component index+fromVectorSv :: Int -- ^ vector dimension+             -> VU.Vector (Int, a) -- ^ vector components (in increasing order)+             -> SVector a+fromVectorSv = SV  -- | 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)]+instance NFData (DVector a)  fromListDv :: VU.Unbox a => [a] -> DVector a fromListDv ll = DV $ VU.fromList ll+fromVectorDv :: VU.Vector a -> DVector a+fromVectorDv = DV toListDv :: (VU.Unbox a) => DVector a -> [a] toListDv (DV v) = VU.toList v @@ -108,11 +133,11 @@ -- 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)+  _rpThreshold ::  !d+  , _rpMargin :: {-# UNPACK #-} !(Margin d)   , _rpL :: !(RPT d a)   , _rpR :: !(RPT d a) }-  | Tip { _rpData :: a }+  | Tip { _rpData :: !a }   deriving (Eq, Show, Generic, Functor, Foldable, Traversable) instance (Serialise a, Serialise d) => Serialise (RPT d a) makeLensesFor [("_rpData", "rpData")] ''RPT@@ -127,14 +152,31 @@ -- 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+  , _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) +-- | A random projection forest is an ordered set of 'RPTree's+--+-- This supports efficient updates of the ensemble in the streaming/online setting. type RPForest d a = IM.IntMap (RPTree d a) +-- | Serialise each tree in the 'RPForest' as a separate bytestring+serialiseRPForest :: (Serialise d, Serialise a, VU.Unbox d) =>+                     RPForest d a+                  -> [LBS.ByteString] -- ^ the order is undefined+serialiseRPForest tt = serialise `map` toList tt++-- | Deserialise each tree in the 'RPForest' from a separate bytestring and reconstruct+deserialiseRPForest :: (Serialise d, Serialise a, VU.Unbox d) =>+                       [LBS.ByteString]+                    -> Either String (RPForest d a) -- ^ the order is undefined+deserialiseRPForest bss = case deserialiseOrFail `traverse` bss of+  Left e -> Left (show e)+  Right xs -> Right $ IM.fromList $ zip [0 ..] xs+ rpTreeData :: Traversal' (RPTree d a) a rpTreeData = rpTree . rpData @@ -152,6 +194,7 @@ -- -- points in 2d -- data P a = P !a !a deriving (Eq, Show) +-- | Scale a vector class Scale v where   (.*) :: (VU.Unbox a, Num a) => a -> v a -> v a instance Scale SVector where@@ -271,7 +314,7 @@ -- | 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 +sumSS = binSS (+) 0  -- | Vector difference diffSD :: (VG.Vector u (Int, a), VG.Vector v a, VU.Unbox a, Num a) =>@@ -302,7 +345,7 @@             GT -> Just ((ir, f z  xr), (i1     , succ i2))  -+-- FIXME the return type of a sparse-dense binary operation depends on the operator itself (S * D = S , S + D = D ), so 'binSD' must be changed 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@@ -318,6 +361,7 @@             y = f xl xr  +{-# SCC partitionAtMedian #-} -- | 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@@ -329,9 +373,7 @@     -- (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+    thr = inns VG.! nh -- inner product threshold,  mgl < thr < mgr     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

@@ -0,0 +1,50 @@+{-# options_ghc -Wno-unused-imports #-}+module Data.RPTree.Internal.Testing where++import Control.Monad.IO.Class (MonadIO(..))+import GHC.Word (Word8, Word64)++-- conduit+import qualified Data.Conduit as C (ConduitT, runConduit, yield, await, transPipe)+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)+-- splitmix+import System.Random.SplitMix (initSMGen, unseedSMGen)+-- splitmix-distributions+import System.Random.SplitMix.Distributions (GenT)++import Data.RPTree.Internal (SVector, fromListSv, DVector, fromListDv)+import Data.RPTree.Gen (dense, sparse, normal2, normalSparse2, normalDense2)+import Data.RPTree.Conduit (dataSource)++data BenchConfig = BenchConfig {+  bcDescription :: String+  , bcMaxTreeDepth :: Int+  , bcMinLeafSize :: Int+  , bcNumTrees :: Int+  , bcChunkSize :: Int+  , bcNZDensity :: Double+  , bcVectorDim :: Int+  , bcDataSize :: Int+  , bcNumQueryPoints :: Int+                               } deriving (Show)++randSeed :: MonadIO m => m Word64+randSeed = liftIO (fst . unseedSMGen <$> initSMGen)+++-- | binary mixture of isotropic Gaussian rvs+datD :: Monad m =>+        Int -- ^ number of data points+     -> Int -- ^ vector dimension+     -> C.ConduitT i (DVector Double) (GenT m) ()+datD n d = dataSource n $ normalDense2 d++-- | binary mixture of isotropic Gaussian rvs with sparse components+datS :: Monad m =>+        Int -- ^ number of data points+     -> Int -- ^ vector dimension+     -> Double -- ^ nonzero density+     -> C.ConduitT i (SVector Double) (GenT m) ()+datS n d pnz = dataSource n $ normalSparse2 pnz d