diff --git a/Changelog.md b/Changelog.md new file mode 100644 --- /dev/null +++ b/Changelog.md @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -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 + diff --git a/bench/time/Main.hs b/bench/time/Main.hs new file mode 100644 --- /dev/null +++ b/bench/time/Main.hs @@ -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)] diff --git a/rp-tree.cabal b/rp-tree.cabal --- a/rp-tree.cabal +++ b/rp-tree.cabal @@ -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 diff --git a/src/Data/RPTree.hs b/src/Data/RPTree.hs --- a/src/Data/RPTree.hs +++ b/src/Data/RPTree.hs @@ -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 diff --git a/src/Data/RPTree/Conduit.hs b/src/Data/RPTree/Conduit.hs --- a/src/Data/RPTree/Conduit.hs +++ b/src/Data/RPTree/Conduit.hs @@ -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 diff --git a/src/Data/RPTree/Gen.hs b/src/Data/RPTree/Gen.hs --- a/src/Data/RPTree/Gen.hs +++ b/src/Data/RPTree/Gen.hs @@ -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 diff --git a/src/Data/RPTree/Internal.hs b/src/Data/RPTree/Internal.hs --- a/src/Data/RPTree/Internal.hs +++ b/src/Data/RPTree/Internal.hs @@ -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 diff --git a/src/Data/RPTree/Internal/Testing.hs b/src/Data/RPTree/Internal/Testing.hs new file mode 100644 --- /dev/null +++ b/src/Data/RPTree/Internal/Testing.hs @@ -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