{-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE ConstrainedClassMethods #-} {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE MultiWayIf #-} {-# LANGUAGE NamedFieldPuns #-} {-# LANGUAGE RecordWildCards #-} module Extrapolation (Extrapolator(..), LinearExtrapolator(..), linearDelay, distanceAlongLine, euclid) where import Data.Foldable (maximumBy, minimumBy) import Data.Function (on) import Data.List.NonEmpty (NonEmpty) import qualified Data.List.NonEmpty as NE import qualified Data.Map as M import Data.Time (Day, UTCTime (UTCTime, utctDay), diffUTCTime, getCurrentTime, nominalDiffTimeToSeconds) import qualified Data.Vector as V import GHC.Float (int2Double) import GHC.IO (unsafePerformIO) import Conduit (MonadIO (liftIO)) import Data.List (sortBy, sortOn) import Data.Ord (Down (..)) import GTFS (Depth (Deep), GTFS (..), Seconds (..), Shape (..), Station (stationName), Stop (..), Time, Trip (..), seconds2Double, stationGeopos, toSeconds) import Persist (Running (..), TrainAnchor (..), TrainPing (..)) import Server.Util (utcToSeconds) -- | Determines how to extrapolate delays (and potentially other things) from the real-time -- data sent in by the OBU. Potentially useful to swap out the algorithm, or give it options. -- TODO: maybe split into two classes? class Extrapolator a where -- | here's a position ping, guess things from that! extrapolateAnchorFromPing :: a -> GTFS -> Running -> TrainPing -> TrainAnchor -- | extrapolate status at some time (i.e. "how much delay does the train have *now*?") extrapolateAtSeconds :: a -> NonEmpty TrainAnchor -> Seconds -> Maybe TrainAnchor -- | extrapolate status at some places (i.e. "how much delay will it have at the next station?") extrapolateAtPosition :: a -> NonEmpty TrainAnchor -> Double -> Maybe TrainAnchor data LinearExtrapolator = LinearExtrapolator instance Extrapolator LinearExtrapolator where extrapolateAtSeconds _ history secondsNow = fmap (minimumBy (compare `on` difference)) $ NE.nonEmpty $ NE.filter (\a -> trainAnchorWhen a < secondsNow) history where difference status = secondsNow - (trainAnchorWhen status) -- note that this sorts (descending) for time first as a tie-breaker -- (in case the train just stands still for a while, take the most recent update) extrapolateAtPosition _ history positionNow = fmap (minimumBy (compare `on` difference)) $ NE.nonEmpty $ sortOn (Down . trainAnchorWhen) $ NE.filter (\a -> trainAnchorSequence a < positionNow) history where difference status = positionNow - (trainAnchorSequence status) extrapolateAnchorFromPing _ gtfs@GTFS{..} Running{..} ping@TrainPing{..} = TrainAnchor { trainAnchorCreated = trainPingTimestamp , trainAnchorTrip = runningTrip , trainAnchorDay = runningDay , trainAnchorWhen = utcToSeconds trainPingTimestamp runningDay , trainAnchorSequence , trainAnchorDelay , trainAnchorMsg = Nothing } where Just trip = M.lookup runningTrip trips (trainAnchorDelay, trainAnchorSequence) = linearDelay gtfs trip ping runningDay linearDelay :: GTFS -> Trip Deep Deep -> TrainPing -> Day -> (Seconds, Double) linearDelay GTFS{..} trip@Trip{..} TrainPing{..} runningDay = (observedDelay, observedSequence) where -- | at which sequence number is the ping? observedSequence = int2Double (stopSequence lastStop) + observedProgress * int2Double (stopSequence nextStop - stopSequence lastStop) -- | how much later/earlier is the ping than would be expected? observedDelay = Seconds $ round $ (expectedProgress - observedProgress) * int2Double (unSeconds expectedTravelTime) -- if the hypothetical on-time train is already at (or past) the next station, -- just add the time distance we're behind + if expectedProgress /= 1 then 0 else seconds2Double (utcToSeconds trainPingTimestamp runningDay - toSeconds (stopArrival nextStop) tzseries runningDay) -- | how far along towards the next station is the ping (between 0 and 1)? observedProgress = distanceAlongLine line (stationGeopos $ stopStation lastStop) closestPoint / distanceAlongLine line (stationGeopos $ stopStation lastStop) (stationGeopos $ stopStation nextStop) -- | to compare: where would a linearly-moving train be (between 0 and 1)? expectedProgress = if | p < 0 -> 0 | p > 1 -> 1 | otherwise -> p where p = seconds2Double (utcToSeconds trainPingTimestamp runningDay - toSeconds (stopDeparture lastStop) tzseries runningDay) / seconds2Double expectedTravelTime -- | how long do we expect the trip from last to next station to take? expectedTravelTime = toSeconds (stopArrival nextStop) tzseries runningDay - toSeconds (stopDeparture lastStop) tzseries runningDay closestPoint = minimumBy (compare `on` euclid (trainPingLat, trainPingLong)) line line = shapePoints tripShape lastStop = tripStops V.! (nextIndex - 1) nextStop = tripStops V.! nextIndex -- | index of the /next/ stop in the list, except when we're already at the last stop -- (in which case it stays the same) nextIndex = if | null remaining -> length tripStops - 1 | idx' == 0 -> 1 | otherwise -> idx' where idx' = fst $ V.minimumBy (compare `on` snd) remaining remaining = V.filter (\(_,dist) -> dist > 0) $ V.indexed $ fmap (distanceAlongLine line closestPoint . stationGeopos . stopStation) tripStops distanceAlongLine :: V.Vector (Double, Double) -> (Double, Double) -> (Double, Double) -> Double distanceAlongLine line p1 p2 = along2 - along1 where along1 = along p1 along2 = along p2 along p@(x,y) = sumSegments $ V.take (index + 1) line where index = V.minIndexBy (compare `on` euclid p) line sumSegments :: V.Vector (Double, Double) -> Double sumSegments line = snd $ foldl (\(p,a) p' -> (p', a + euclid p p')) (V.head line,0) $ line -- | euclidean distance. Notably not applicable when you're on a sphere -- (but good enough when the sphere is the earth) euclid :: Floating f => (f,f) -> (f,f) -> f euclid (x1,y1) (x2,y2) = sqrt (x*x + y*y) where x = x1 - x2 y = y1 - y2