{-# LANGUAGE TupleSections #-}
module Data.GPS.Trail
( -- * Types
AvgMethod(..)
, Selected(..)
, PointGrouping
, TransformGrouping
-- * Utility Functions
, isSelected
, isNotSelected
, onSelected
, selLength
-- * Trail Functions
-- ** Queries
, totalDistance
, totalTime
, avgSpeeds
, slidingAverageSpeed
, closestDistance
, convexHull
-- ** Transformations
, bezierCurveAt
, bezierCurve
, linearTime
, filterPoints
-- ** Grouping Methods
, betweenSpeeds
, restLocations
, spansTime
, everyNPoints
-- ** Group Transformations
, intersectionOf
, invertSelection
, firstGrouping
, lastGrouping
, unionOf
, refineGrouping
, (/\), (\/)
-- ** Composite Operations (Higher Level)
, smoothRests
, smoothSome
, smoothMore
-- * Misc
, bezierPoint
) where
import Text.Show.Functions ()
import Data.GPS.Core hiding (fix)
import Text.XML.XSD.DateTime (fromUTCTime)
import Control.Arrow (first, second)
import Control.Monad
import Data.Fixed (mod')
import Data.Function (on,fix)
import Data.List as L
import Data.Maybe
import Data.Ord
import Data.Time
import Data.Lens.Common
import Statistics.Function as F
import Statistics.Sample
import qualified Data.Vector.Unboxed as V
takeWhileEnd :: (a -> Bool) -> [a] -> (Maybe a, [a],[a])
takeWhileEnd p xs = go xs Nothing
where
go [] e = (e, [], [])
go (a:as) e
| p a = let (e',xs,zs) = go as (Just a) in (e',a:xs,zs)
| otherwise = (e,[], a:as)
data AvgMethod c
= AvgMean -- ^ Obtain the 'mean' of the considered points
| AvgHarmonicMean -- ^ Obtain the 'harmonicMean'
| AvgGeometricMean -- ^ Obtain the 'geometricMean'
| AvgMedian -- ^ Obtain the median of the considered points
| AvgEndPoints -- ^ Compute the speed considering only the given endpoints
| AvgMinOf [AvgMethod c] -- ^ Take the minimum of the speeds from the given methods
| AvgWith ([c] -> Speed)
-- | @avgSpeeds n points@
-- Average speed using a window of up to @n@ seconds and averaging by taking the
-- Median ('AvgMedian').
avgSpeeds :: (Coordinate a, TimeL a) => NominalDiffTime -> Trail a -> [(UTCTime, Speed)]
avgSpeeds = slidingAverageSpeed AvgHarmonicMean
-- | @slidingAverageSpeed m n@ Average speed using a moving window of up to @n@ seconds
-- and an 'AvgMethod' of @m@.
slidingAverageSpeed :: (Coordinate a, TimeL a) =>
AvgMethod a -> NominalDiffTime -> Trail a -> [(UTCTime, Speed)]
slidingAverageSpeed _ _ [] = []
slidingAverageSpeed m minTime xs =
let pts = map unSelect (spansTime minTime xs)
spds = map (getAvg m) pts
times = map getAvgTimes pts
in concatMap maybeToList $ zipWith (\t s -> fmap (,s) t) times spds
where
getTimeDiff a b = on (liftM2 diffUTCTime) getUTCTime a b
-- getAvg :: [] -> AvgMethod -> Speed
getAvg _ [] = 0
getAvg _ [x] = 0
getAvg m cs =
let ss = getSpeedsV cs
in case m of
AvgMean -> mean ss
AvgHarmonicMean -> harmonicMean ss
AvgGeometricMean -> geometricMean ss
AvgMedian ->
let ss' = F.sort $ getSpeedsV cs
len = V.length ss'
mid = len `div` 2
in if V.length ss' < 3
then mean ss'
else if odd len then ss' V.! mid else mean (V.slice mid 2 ss')
AvgEndPoints -> fromMaybe 0 $ speed (head cs) (last cs)
AvgMinOf as -> minimum $ map (flip getAvg cs) as
AvgWith f -> f cs
getAvgTimes [] = Nothing
getAvgTimes [x] = getUTCTime x
getAvgTimes ps = getAvgTime (head ps) (last ps)
getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (getUTCTime a)
getSpeedsV = V.fromList . getSpeeds
getSpeeds zs = concatMap maybeToList $ zipWith speed zs (drop 1 zs)
-- | A PointGrouping is a function that selects segments of a trail.
--
-- Grouping point _does not_ result in deleted points. It is always true that:
--
-- forall g :: PointGrouping c -->
-- concatMap unSelect (g ts) == ts
--
-- The purpose of grouping is usually for later processing. Any desire to drop
-- points that didn't meet a particular grouping criteria can be filled with
-- a composition with 'filter' (or directly via 'filterPoints').
type PointGrouping c = Trail c -> [Selected (Trail c)]
-- | Given a selection of coordinates, transform the selected
-- coordinates in some way (while leaving the non-selected
-- coordinates unaffected).
type TransformGrouping c = [Selected (Trail c)] -> [Selected (Trail c)]
-- | When grouping points, lists of points are either marked as 'Select' or 'NotSelect'.
data Selected a = Select {unSelect :: a} | NotSelect {unSelect :: a}
deriving (Eq, Ord, Show)
isSelected :: Selected a -> Bool
isSelected (Select _) = True
isSelected _ = False
isNotSelected :: Selected a -> Bool
isNotSelected = not . isSelected
selLength :: Selected [a] -> Int
selLength = length . unSelect
onSelected :: (a -> b) -> (a -> b) -> Selected a -> b
onSelected f _ (Select a) = f a
onSelected _ g (NotSelect a) = g a
instance Functor Selected where
fmap f (Select x) = Select $ f x
fmap f (NotSelect x) = NotSelect $ f x
dropExact :: Int -> [Selected [a]] -> [Selected [a]]
dropExact i [] = []
dropExact i (x:xs) =
case compare (selLength x) i of
EQ -> xs
LT -> dropExact (i - selLength x) xs
GT -> fmap (drop i) x : xs
-- | Groups trail segments into contiguous points within the speed
-- and all others outside of the speed. The "speed" from point p(i)
-- to p(i+1) is associated with p(i) (execpt for the first speed
-- value, which is associated with both the first and second point)
betweenSpeeds :: (Coordinate a, TimeL a) => Double -> Double -> PointGrouping a
betweenSpeeds low hi ps =
let spds = concatMap maybeToList $ zipWith speed ps (drop 1 ps)
psSpds = [(p,s) | p <- ps, s <- maybeToList (listToMaybe spds) ++ spds]
inRange x = x >= low && x <= hi
chunk [] = []
chunk xs@(x:_) =
let op p = if inRange (snd x) then first Select . span p else first NotSelect . break p
(r,rest) = op (inRange . snd) xs
in r : chunk xs
in map (fmap (map fst)) $ chunk psSpds
-- | A "rest point" means the coordinates remain within a given distance
-- for at least a particular amount of time.
restLocations :: (Coordinate a, TimeL a) => Distance -> NominalDiffTime -> PointGrouping a
restLocations d s ps =
let consToFirst x [] = [NotSelect [x]]
consToFirst x (a:as) = (fmap (x:) a) : as
go [] [] = []
go [] nonRests = [NotSelect $ reverse nonRests]
go (a:as) nonRests =
case takeWhileEnd ((<=) d . distance a) as of
(Just l, close, far) ->
case (getUTCTime a, getUTCTime l) of
(Just t1, Just t2) ->
let diff = diffUTCTime t2 t1
in if diff >= s then NotSelect (reverse nonRests) : Select (a:close) : go far [] else go as (a:nonRests)
_ -> consToFirst a $ go as nonRests
_ -> consToFirst a $ go as nonRests
in go ps []
-- |chunking points into groups spanning at most the given time
-- interval.
spansTime :: (Coordinate a, TimeL a) => NominalDiffTime -> PointGrouping a
spansTime n ps =
let times = mkTimePair ps
chunk [] = []
chunk xs@(x:_) =
let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs
in if null good then [xs] else good : chunk rest
in map (Select . map fst) $ chunk times
-- | intersects the given groupings
intersectionOf :: (Coordinate a, TimeL a) => [PointGrouping a] -> PointGrouping a
intersectionOf gs ps =
let groupings = map ($ ps) gs
-- chunk :: [[Selected [pnts]]] -> pnts -> [pnts]
chunk _ [] = []
chunk ggs xs =
let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs -- FIXME this is all manner of broken
sel = if all isSelected (concatMap (take 1) ggs) then Select else NotSelect
(c,rest) = splitAt minLen xs
in sel c : chunk (filter (not . null) $ map (dropExact minLen) ggs) rest
in chunk groupings ps
-- | Union all the groupings
unionOf :: (Coordinate a, TimeL a) => [PointGrouping a] -> PointGrouping a
unionOf gs ps =
let groupings = map ($ ps) gs
chunk _ [] = []
chunk ggs xs =
let getSegs = concatMap (take 1)
segs = getSegs ggs
len =
if any isSelected segs
then max 1 . maximum . getSegs . map (map selLength) . map (filter isSelected) $ ggs
else max 1 . minimum . getSegs . map (map selLength) $ ggs
sel = if any isSelected segs then Select else NotSelect
(c,rest) = splitAt len xs
in sel c : chunk (filter (not . null) $ map (dropExact len) ggs) rest
in chunk groupings ps
-- | Intersection binary operator
(/\) :: [Selected (Trail a)] -> TransformGrouping a
(/\) xs [] = []
(/\) [] ys = []
(/\) xsL@(Select x:xs) ysL@(Select y:ys) =
let z = if length x < length y then x else y
xs' = selListDrop (length z) xsL
ys' = selListDrop (length z) ysL
in Select z : (xs' /\ ys')
(/\) xs (NotSelect y:ys) = NotSelect y : (selListDrop (length y) xs /\ ys)
(/\) (NotSelect x:xs) ys = NotSelect x : (xs /\ selListDrop (length x) ys)
-- | Union binary operator
(\/) :: [Selected (Trail a)] -> TransformGrouping a
(\/) xs [] = xs
(\/) [] ys = ys
(\/) (Select x:xs) (Select y : ys) =
let xLen = length x
yLen = length y
in if xLen < yLen
then (Select y :) (selListDrop (yLen - xLen) xs \/ ys)
else (Select x :) (xs \/ selListDrop (xLen - yLen) ys)
(\/) (Select x:xs) ys = Select x : selListDrop (length x) ys
(\/) xs (Select y:ys) = Select y : selListDrop (length y) xs
(\/) xsL@(NotSelect x:xs) ysL@(NotSelect y:ys) =
let xLen = length x
yLen = length y
in if xLen < yLen
then (NotSelect x:) (xs \/ selListDrop xLen ysL)
else (NotSelect y:) (selListDrop yLen xsL \/ ys)
selListTake :: Int -> [Selected [a]] -> [Selected [a]]
selListTake 0 _ = []
selListTake n [] = []
selListTake n (x:xs) =
let x' = take n (unSelect x)
in fmap (const x') x : selListTake (n - length x') xs
selListDrop :: Int -> [Selected [a]] -> [Selected [a]]
selListDrop 0 xs = xs
selListDrop n [] = []
selListDrop n (x:xs) =
let x' = drop n (unSelect x)
in fmap (const x') x : selListDrop (n - (selLength x - length x')) xs
-- |Inverts the selected/nonselected segments
invertSelection :: TransformGrouping a
invertSelection = map (onSelected NotSelect Select)
-- |@firstGrouping f ps@ only the first segment remains 'Select'ed, and only
-- if it was already selected by @f@.
firstGrouping :: TransformGrouping a
firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps)
-- | Only the last segment, if any, is selected (note: the current
-- implementation is inefficient, using 'reverse')
lastGrouping :: TransformGrouping a
lastGrouping ps = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps')
-- | chunk the trail into groups of N points
everyNPoints :: Int -> PointGrouping a
everyNPoints n ps
| n <= 0 = [NotSelect ps]
| otherwise = go ps
where
go [] = []
go xs = let (h,t) = splitAt n xs in Select h : go t
-- |For every selected group, refine the selection using the second
-- grouping method. This differs from 'IntersectionOf' by restarting
-- the second grouping algorithm at the beginning each group selected
-- by the first algorithm.
refineGrouping :: PointGrouping a -> TransformGrouping a
refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x]))
-- |Remove all points that remain 'NotSelect'ed by the given grouping algorithm.
filterPoints :: PointGrouping a -> Trail a -> Trail a
filterPoints g = concatMap unSelect . filter isSelected . g
-- Extract the time from each coordinate. If no time is available then
-- the coordinate is dropped!
mkTimePair :: (TimeL a) => Trail a -> [(a,UTCTime)]
mkTimePair xs =
let timesM = map (\x-> fmap (x,) $ getUTCTime x) xs
in concatMap maybeToList timesM
-- |Construct a bezier curve using the provided trail. Construct a
-- new trail by sampling the given bezier curve at the given times.
-- The current implementation assumes the times of the input
-- coordinates are available and all equal (Ex: all points are 5
-- seconds apart), the results will be poor if this is not the case!
bezierCurveAt :: (Coordinate a, TimeL a) => [UTCTime] -> Trail a -> Trail a
bezierCurveAt _ [] = []
bezierCurveAt selectedTimes xs =
let timesDef = mkTimePair xs
end = last timesDef
top = head timesDef
totalTime = diffUTCTime (snd end) (snd top)
times = if null selectedTimes then map snd timesDef else selectedTimes
diffTimes = [diffUTCTime t (snd top) / totalTime | t <- times]
queryTimes = map realToFrac diffTimes
in if totalTime <= 0 || any (\x -> x < 0 || x > 1) queryTimes
then xs -- error "bezierCurveAt has a out-of-bound time!"
else
if null timesDef || any (\x -> x < 0 || x > 1) queryTimes
then xs
else let curvePoints = (map (bezierPoint xs) queryTimes)
newTimes = [addUTCTime t (snd top) | t <- diffTimes]
in zipWith ((timeL ^=) . Just . fromUTCTime) newTimes curvePoints
bezierPoint :: (Coordinate a) => [a] -> Double -> a
bezierPoint pnts t = go pnts
where
go [] = error "GPS Package: Can not create a bezier point from an empty list"
go [p] = p
go ps = interpolate (go (init ps)) (go (tail ps)) t
-- |Interpolate selected points onto a bezier curve. Note this gets
-- exponentially more expensive with the length of the segement being
-- transformed - it is not advisable to perform this operation on
-- trail segements with more than ten points!
bezierCurve :: (Coordinate a, TimeL a) => [Selected (Trail a)] -> Trail a
bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id)
-- |Filters out any points that go backward in time (thus must not be
-- valid if this is a trail)
linearTime :: (Coordinate a, TimeL a) => [a] -> [a]
linearTime [] = []
linearTime (p:ps) = go (getUTCTime p) ps
where
go _ [] = []
go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps
-- |Returns the closest distance between two trails (or Nothing if a
-- trail is empty). Inefficient implementation:
-- O( (n * m) * log (n * m) )
closestDistance :: (Coordinate a) => Trail a -> Trail a -> Maybe Distance
closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]
-- | Find the total distance traveled
totalDistance :: (Coordinate a) => [a] -> Distance
totalDistance as = sum $ zipWith distance as (drop 1 as)
totalTime :: TimeL a => Trail a -> NominalDiffTime
totalTime [] = 0
totalTime xs@(x:_) = fromMaybe 0 $ liftM2 diffUTCTime (getUTCTime x) (getUTCTime $ last xs)
-- | Uses Grahams scan to compute the convex hull of the given points.
-- This operation requires sorting of the points, so don't try it unless
-- you have notably more memory than the list of points will consume.
convexHull :: (Eq c, Coordinate c) => [c] -> [c]
convexHull xs =
let first = southMost xs
in case first of
Nothing -> []
Just f ->
let sorted = L.sortBy (comparing (eastZeroHeading f)) (filter (/= f) xs)
in case sorted of
(a:b:cs) -> grahamScan (b:a:f:[]) cs
cs -> f : cs
where
grahamScan [] _ = []
grahamScan ps [] = ps
grahamScan (x:[]) _ = [x]
grahamScan (p2:p1:ps) (x:xs) =
case turn p1 p2 x of
LeftTurn -> grahamScan (x:p2:p1:ps) xs
Straight -> grahamScan (x:p2:p1:ps) xs
_ -> grahamScan (p1:ps) (x:xs)
eastZeroHeading :: (Coordinate c) => c -> c -> Heading
eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s
data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)
turn :: (Coordinate c) => c -> c -> c -> Turn
turn a b c =
let h1 = eastZeroHeading a b
h2 = eastZeroHeading b c
d = h2 - h1
in if d >= 0 && d < pi then LeftTurn else RightTurn
-- | Find the southmost point
southMost :: (LatL c) => [c] -> Maybe c
southMost [] = Nothing
southMost cs = Just . minimumBy (comparing (^. latL)) $ cs
---------- COMPOSIT OPERATIONS ---------------
-- These operations are simply implemented using the previously
-- defined functions. They can serve either for concise use for novice
-- users or as instructional examples.
------------------------------------------
smoothRests :: (Coordinate a, TimeL a) => Trail a -> Trail a
smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60
smoothSome :: (Coordinate a, TimeL a) => Trail a -> Trail a
smoothSome = gSmoothSome 7
gSmoothSome n = bezierCurve . everyNPoints n
gSmoothMore n k ps =
let ps' = gSmoothSome n ps
(h,t) = splitAt k ps'
in h ++ gSmoothSome n t
smoothMore :: (Coordinate a, TimeL a) => Trail a -> Trail a
smoothMore
= gSmoothMore 7 3 . gSmoothMore 3 1 . gSmoothMore 5 2
. gSmoothMore 3 1 . gSmoothMore 5 2 . gSmoothMore 7 3
slidingWindow :: Int -> Int -> ([a] -> [a]) -> [a] -> [a]
slidingWindow width step f trail = go trail
where
go [] = []
go xs =
let (window,e) = splitAt width xs
(hT,eT) = splitAt step (f window)
in hT ++ go (eT ++ e)