Algorithms.Geometry.WellSeparatedPairDecomposition.WSPD

Documentation

fairSplitTree :: (Fractional r, Ord r, Arity d, 1 <= d, Show r, Show p) => NonEmpty (Point d r :+ p) -> SplitTree d p r ()

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]

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)

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))

Assign the points to their the correct class. The Nothing class is considered the last class

transpose :: Arity d => Vector d (Vector a) -> Vector (Vector d a)

distributePoints'

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

reIndexPoints :: (Arity d, 1 <= d) => Vector d (PointSeq d (Idx :+ p) r) -> Vector d (PointSeq d (Idx :+ p) r)

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)

ST monad with access to the vector storign the level of the points.

assignLevels

Arguments

:: (Fractional r, Ord r, Arity d, KnownNat 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))

Remove already assigned pts from the ends of all vectors.

assignLevel :: (c :+ (Idx :+ p)) -> Level -> RST s ()

Assign level l to point p

levelOf :: (c :+ (Idx :+ p)) -> RST s (Maybe Level)

Get the level of a point

hasLevel :: (c :+ (Idx :+ p)) -> RST s Bool

Test if the point already has a level assigned to it.

compactEnds' :: PointSeq d (Idx :+ p) r -> RST s (PointSeq d (Idx :+ p) r)

Remove allready assigned points from the sequence

pre: there are points remaining

findAndCompact

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

Find the widest dimension of the point set

pre: points are sorted according to their dimension

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)

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]

areWellSeparated

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

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

boxBox1 :: (Floating r, Ord r, Arity d) => r -> Box d p r -> Box d p r -> Bool

Helper stuff

maxWidth :: (Arity d, KnownNat d, Num r) => SplitTree d p r a -> r

Computes the maximum width of a splitTree

bbOf :: Ord r => SplitTree d p r a -> Box d () r

Computes the bounding box of a split tree

children' :: BinLeafTree v a -> [BinLeafTree v a]

fromSeqUnsafe :: Seq a -> ViewL1 a

ix' :: (Arity d, KnownNat d) => Int -> Lens' (Vector d a) a

Turn a traversal into lens

dropIdx :: (core :+ (t :+ extra)) -> core :+ extra