aboutsummaryrefslogtreecommitdiff
path: root/lib/Extrapolation.hs
blob: 759b31ef28d4022bcf1e1c6992d23d85f59af009 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE MultiWayIf          #-}
{-# 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           API                                 (SentPing (..))
import           Conduit                             (MonadIO (liftIO))
import           Data.List                           (sortBy, sortOn)
import           Data.Ord                            (Down (..))
import           Data.Time.LocalTime.TimeZone.Series (TimeZoneSeries)
import           GTFS                                (Seconds (..),
                                                      seconds2Double, toSeconds)
import           Persist                             (Geopos (..),
                                                      ShapePoint (shapePointGeopos),
                                                      Station (..), Stop (..),
                                                      Ticket (..), TicketId,
                                                      Token (..), Tracker (..),
                                                      TrainAnchor (..))
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 strategy where
  -- | here's a position ping, guess things from that!
  extrapolateAnchorFromPing
    :: strategy
    -> TicketId
    -> Ticket
    -> V.Vector (Stop, Station, TimeZoneSeries)
    -> V.Vector ShapePoint
    -> SentPing
    -> TrainAnchor

  -- | extrapolate status at some time (i.e. "how much delay does the train have *now*?")
  extrapolateAtSeconds :: strategy -> NonEmpty TrainAnchor -> Seconds -> Maybe TrainAnchor
  -- | extrapolate status at some places (i.e. "how much delay will it have at the next station?")
  extrapolateAtPosition :: strategy -> 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 . trainAnchorCreated)
    $ NE.filter (\a -> trainAnchorSequence a < positionNow) history
    where difference status = positionNow - trainAnchorSequence status

  extrapolateAnchorFromPing _ ticketId Ticket{..} stops shape ping@SentPing{..} = TrainAnchor
    { trainAnchorCreated = sentPingTimestamp
    , trainAnchorTicket = ticketId
    , trainAnchorWhen = utcToSeconds sentPingTimestamp ticketDay
    , trainAnchorSequence
    , trainAnchorDelay
    , trainAnchorMsg = Nothing
    }
    where
          (trainAnchorDelay, trainAnchorSequence) = linearDelay stops shape ping ticketDay
          tzseries = undefined

linearDelay :: V.Vector (Stop, Station, TimeZoneSeries) -> V.Vector ShapePoint -> SentPing -> Day -> (Seconds, Double)
linearDelay tripStops shape SentPing{..} runningDay = (observedDelay, observedSequence)
  where -- at which (fractional) 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 expectedProgress == 1
          -- if the hypothetical on-time train is already at (or past) the next station,
          -- just add the time distance we're behind
          then seconds2Double (utcToSeconds sentPingTimestamp runningDay - nextSeconds)
          -- otherwise the above is sufficient
          else 0

        -- how far along towards the next station is the ping (between 0 and 1)?
        observedProgress =
            distanceAlongLine line (stationGeopos lastStation) closestPoint
            / distanceAlongLine line (stationGeopos lastStation) (stationGeopos nextStation)

        -- 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 sentPingTimestamp runningDay - lastSeconds)
                  / seconds2Double expectedTravelTime

        -- scheduled duration between last and next stops
        expectedTravelTime = nextSeconds - lastSeconds

        -- closest point on the shape; this is where we assume the train to be
        closestPoint = minimumBy (compare `on` euclid sentPingGeopos) line

        -- scheduled departure at last & arrival at next stop
        lastSeconds = toSeconds (stopDeparture lastStop) lastTzSeries runningDay
        nextSeconds = toSeconds (stopArrival nextStop) nextTzSeries runningDay

        (lastStop, lastStation, lastTzSeries) = tripStops V.! (nextIndex - 1)
        (nextStop, nextStation, nextTzSeries) = tripStops V.! nextIndex

        line = fmap shapePointGeopos shape

        -- 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 . \(_,stop,_) -> stop) tripStops

-- | approximate (but euclidean) distance along a geoline
distanceAlongLine :: V.Vector Geopos -> Geopos -> Geopos -> Double
distanceAlongLine line p1 p2 = along2 - along1
  where along1 = along p1
        along2 = along p2
        along p@(Geopos (x,y)) =
          sumSegments
          $ V.take (index + 1) line
          where index = V.minIndexBy (compare `on` euclid p) line
        sumSegments :: V.Vector Geopos -> 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 :: Geopos -> Geopos -> Double
euclid (Geopos (x1,y1)) (Geopos (x2,y2)) = sqrt (x*x + y*y)
  where x = x1 - x2
        y = y1 - y2