rp-tree (empty) → 0.1.0.0
raw patch · 11 files changed
+1427/−0 lines, 11 filesdep +QuickCheckdep +basedep +boxessetup-changed
Dependencies added: QuickCheck, base, boxes, bytestring, conduit, containers, deepseq, exceptions, hspec, microlens, microlens-th, mtl, rp-tree, serialise, splitmix-distributions, transformers, vector, vector-algorithms
Files
- LICENSE +30/−0
- README.md +5/−0
- Setup.hs +2/−0
- app/Main.hs +123/−0
- rp-tree.cabal +77/−0
- src/Data/RPTree.hs +359/−0
- src/Data/RPTree/Conduit.hs +215/−0
- src/Data/RPTree/Draw.hs +86/−0
- src/Data/RPTree/Gen.hs +149/−0
- src/Data/RPTree/Internal.hs +380/−0
- test/Spec.hs +1/−0
+ LICENSE view
@@ -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.
+ README.md view
@@ -0,0 +1,5 @@+# rp-tree++++Random projection trees
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/Main.hs view
@@ -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
+ rp-tree.cabal view
@@ -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++
+ src/Data/RPTree.hs view
@@ -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
+ src/Data/RPTree/Conduit.hs view
@@ -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 -- .|++
+ src/Data/RPTree/Draw.hs view
@@ -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]
+ src/Data/RPTree/Gen.hs view
@@ -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)
+ src/Data/RPTree/Internal.hs view
@@ -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
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}