| Copyright | (C) Frank Staals |
|---|---|
| License | see the LICENSE file |
| Maintainer | Frank Staals |
| Safe Haskell | None |
| Language | Haskell2010 |
Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD
Contents
Description
Algorithm to construct a well separated pair decomposition (wspd).
Synopsis
- fairSplitTree :: (Fractional r, Ord r, Arity d, 1 <= d, Show r, Show p) => NonEmpty (Point d r :+ p) -> SplitTree d p r ()
- wellSeparatedPairs :: (Floating r, Ord r, Arity d, Arity (d + 1)) => r -> SplitTree d p r a -> [WSP d p r a]
- fairSplitTree' :: (Fractional r, Ord r, Arity d, 1 <= d, Show r, Show p) => Int -> Vector d (PointSeq d (Idx :+ p) r) -> BinLeafTree Int (Point d r :+ p)
- distributePoints :: (Arity d, Show r, Show p) => Int -> Vector (Maybe Level) -> Vector d (PointSeq d (Idx :+ p) r) -> Vector (Vector d (PointSeq d (Idx :+ p) r))
- transpose :: Arity d => Vector d (Vector a) -> Vector (Vector d a)
- distributePoints' :: Int -> Vector (Maybe Level) -> PointSeq d (Idx :+ p) r -> Vector (PointSeq d (Idx :+ p) r)
- fromSeqUnsafe :: Seq a -> LSeq n a
- reIndexPoints :: (Arity d, 1 <= d) => Vector d (PointSeq d (Idx :+ p) r) -> Vector d (PointSeq d (Idx :+ p) r)
- type RST s = ReaderT (MVector s (Maybe Level)) (ST s)
- assignLevels :: (Fractional r, Ord r, Arity d, Show r, Show p) => Int -> Int -> Vector d (PointSeq d (Idx :+ p) r) -> Level -> [Level] -> RST s (NonEmpty Level)
- compactEnds :: Arity d => Vector d (PointSeq d (Idx :+ p) r) -> RST s (Vector d (PointSeq d (Idx :+ p) r))
- assignLevel :: (c :+ (Idx :+ p)) -> Level -> RST s ()
- levelOf :: (c :+ (Idx :+ p)) -> RST s (Maybe Level)
- hasLevel :: (c :+ (Idx :+ p)) -> RST s Bool
- compactEnds' :: PointSeq d (Idx :+ p) r -> RST s (PointSeq d (Idx :+ p) r)
- findAndCompact :: (Ord r, Arity d, Show r, Show p) => Int -> PointSeq d (Idx :+ p) r -> r -> RST s (PointSeq d (Idx :+ p) r, PointSeq d (Idx :+ p) r)
- widestDimension :: (Num r, Ord r, Arity d) => Vector d (PointSeq d p r) -> Int
- widths :: (Num r, Arity d) => Vector d (PointSeq d p r) -> Vector d r
- extends :: Arity d => Vector d (PointSeq d p r) -> Vector d (Range r)
- findPairs :: (Floating r, Ord r, Arity d, Arity (d + 1)) => r -> SplitTree d p r a -> SplitTree d p r a -> [WSP d p r a]
- areWellSeparated :: (Arity d, Arity (d + 1), Fractional r, Ord r) => r -> SplitTree d p r a -> SplitTree d p r a -> Bool
- boxBox :: (Fractional r, Ord r, Arity d, Arity (d + 1)) => r -> Box d p r -> Box d p r -> Bool
- areWellSeparated' :: (Floating r, Ord r, Arity d) => r -> SplitTree d p r a -> SplitTree d p r a -> Bool
- boxBox1 :: (Floating r, Ord r, Arity d) => r -> Box d p r -> Box d p r -> Bool
- maxWidth :: (Arity d, Num r) => SplitTree d p r a -> r
- bbOf :: Ord r => SplitTree d p r a -> Box d () r
- children' :: BinLeafTree v a -> [BinLeafTree v a]
- ix' :: (Arity d, KnownNat d) => Int -> Lens' (Vector d a) a
- dropIdx :: (core :+ (t :+ extra)) -> core :+ extra
Documentation
fairSplitTree :: (Fractional r, Ord r, Arity d, 1 <= d, Show r, Show p) => NonEmpty (Point d r :+ p) -> SplitTree d p r () Source #
Construct a split tree
running time: \(O(n \log n)\)
wellSeparatedPairs :: (Floating r, Ord r, Arity d, Arity (d + 1)) => r -> SplitTree d p r a -> [WSP d p r a] Source #
Given a split tree, generate the Well separated pairs
running time: \(O(s^d n)\)
Building the split tree
fairSplitTree' :: (Fractional r, Ord r, Arity d, 1 <= d, Show r, Show p) => Int -> Vector d (PointSeq d (Idx :+ p) r) -> BinLeafTree Int (Point d r :+ p) Source #
Given the points, sorted in every dimension, recursively build a split tree
The algorithm works in rounds. Each round takes O(n) time, and halves the number of points. Thus, the total running time is O(n log n).
The algorithm essentially builds a path in the split tree; at every node on the path that we construct, we split the point set into two sets (L,R) according to the longest side of the bounding box.
The smaller set is "assigned" to the current node and set asside. We continue to build the path with the larger set until the total number of items remaining is less than n/2.
To start the next round, each node on the path needs to have the points assigned to that node, sorted in each dimension (i.e. the Vector (PointSeq))'s. Since we have the level assignment, we can compute these lists by traversing each original input list (i.e. one for every dimension) once, and partition the points based on their level assignment.
distributePoints :: (Arity d, Show r, Show p) => Int -> Vector (Maybe Level) -> Vector d (PointSeq d (Idx :+ p) r) -> Vector (Vector d (PointSeq d (Idx :+ p) r)) Source #
Assign the points to their the correct class. The Nothing class is
considered the last class
Arguments
| :: Int | number of classes |
| -> Vector (Maybe Level) | level assignment |
| -> PointSeq d (Idx :+ p) r | input points |
| -> Vector (PointSeq d (Idx :+ p) r) |
Assign the points to their the correct class. The Nothing class is
considered the last class
fromSeqUnsafe :: Seq a -> LSeq n a Source #
reIndexPoints :: (Arity d, 1 <= d) => Vector d (PointSeq d (Idx :+ p) r) -> Vector d (PointSeq d (Idx :+ p) r) Source #
Given a sequence of points, whose index is increasing in the first dimension, i.e. if idx p < idx q, then p[0] < q[0]. Reindex the points so that they again have an index in the range [0,..,n'], where n' is the new number of points.
running time: O(n' * d) (more or less; we are actually using an intmap for the lookups)
alternatively: I can unsafe freeze and thaw an existing vector to pass it
along to use as mapping. Except then I would have to force the evaluation
order, i.e. we cannot be in reIndexPoints for two of the nodes at the same
time.
so, basically, run reIndex points in ST as well.
type RST s = ReaderT (MVector s (Maybe Level)) (ST s) Source #
ST monad with access to the vector storign the level of the points.
Arguments
| :: (Fractional r, Ord r, Arity d, Show r, Show p) | |
| => Int | Number of items we need to collect |
| -> Int | Number of items we collected so far |
| -> Vector d (PointSeq d (Idx :+ p) r) | |
| -> Level | next level to use |
| -> [Level] | Levels used so far |
| -> RST s (NonEmpty Level) |
Assigns the points to a level. Returns the list of levels used. The first level in the list is the level assigned to the rest of the nodes. Their level is actually still set to Nothing in the underlying array.
compactEnds :: Arity d => Vector d (PointSeq d (Idx :+ p) r) -> RST s (Vector d (PointSeq d (Idx :+ p) r)) Source #
Remove already assigned pts from the ends of all vectors.
hasLevel :: (c :+ (Idx :+ p)) -> RST s Bool Source #
Test if the point already has a level assigned to it.
compactEnds' :: PointSeq d (Idx :+ p) r -> RST s (PointSeq d (Idx :+ p) r) Source #
Remove allready assigned points from the sequence
pre: there are points remaining
Arguments
| :: (Ord r, Arity d, Show r, Show p) | |
| => Int | the dimension we are in, i.e. so that we know which coordinate of the point to compare |
| -> PointSeq d (Idx :+ p) r | |
| -> r | the mid point |
| -> RST s (PointSeq d (Idx :+ p) r, PointSeq d (Idx :+ p) r) |
Given the points, ordered by their j^th coordinate, split the point set into a "left" and a "right" half, i.e. the points whose j^th coordinate is at most the given mid point m, and the points whose j^th coordinate is larger than m.
We return a pair (Largest set, Smallest set)
fi ndAndCompact works by simultaneously traversing the points from left to right, and from right to left. As soon as we find a point crossing the mid point we stop and return. Thus, in principle this takes only O(|Smallest set|) time.
running time: O(|Smallest set|) + R, where R is the number of *old* points (i.e. points that should have been removed) in the list.
widestDimension :: (Num r, Ord r, Arity d) => Vector d (PointSeq d p r) -> Int Source #
Find the widest dimension of the point set
pre: points are sorted according to their dimension
extends :: Arity d => Vector d (PointSeq d p r) -> Vector d (Range r) Source #
get the extends of the set of points in every dimension, i.e. the left and right boundaries.
pre: points are sorted according to their dimension
Finding Well Separated Pairs
findPairs :: (Floating r, Ord r, Arity d, Arity (d + 1)) => r -> SplitTree d p r a -> SplitTree d p r a -> [WSP d p r a] Source #
Arguments
| :: (Arity d, Arity (d + 1), Fractional r, Ord r) | |
| => r | separation factor |
| -> SplitTree d p r a | |
| -> SplitTree d p r a | |
| -> Bool |
Test if the two sets are well separated with param s
boxBox :: (Fractional r, Ord r, Arity d, Arity (d + 1)) => r -> Box d p r -> Box d p r -> Bool Source #
Test if the two boxes are sufficiently far appart
Alternative def if wellSeparated that uses fractional
areWellSeparated' :: (Floating r, Ord r, Arity d) => r -> SplitTree d p r a -> SplitTree d p r a -> Bool Source #
Helper stuff
maxWidth :: (Arity d, Num r) => SplitTree d p r a -> r Source #
Computes the maximum width of a splitTree
children' :: BinLeafTree v a -> [BinLeafTree v a] Source #