aboutsummaryrefslogtreecommitdiff
path: root/lib/Extrapolation.hs
blob: 770d4ce8a44035fe2a688f929395fa70ae7cd17c (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
{-# LANGUAGE AllowAmbiguousTypes        #-}
{-# LANGUAGE ConstrainedClassMethods    #-}
{-# LANGUAGE ConstraintKinds            #-}
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE LambdaCase                 #-}
{-# LANGUAGE NamedFieldPuns             #-}
{-# LANGUAGE RecordWildCards            #-}

module Extrapolation (Extrapolator(..), LinearExtrapolator, linearDelay, secondsNow) 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)
import           GTFS               (Depth (Deep), GTFS (..), Seconds (..),
                                     Shape (..), Stop (..), Time, Trip (..),
                                     seconds2Double, stationGeopos, toSeconds)
import           Persist            (Running (..), TrainAnchor (..),
                                     TrainPing (..))

class Extrapolator a where
  -- | here's a position ping, guess things from that!
  extrapolateAnchorFromPing :: GTFS -> Running -> TrainPing -> TrainAnchor

  -- | extrapolate status at some time (i.e. "how much delay does the train have *now*?")
  extrapolateAtSeconds :: NonEmpty TrainAnchor -> Seconds -> Maybe TrainAnchor
  -- | extrapolate status at some places (i.e. "how much delay will it have at the next station?")
  extrapolateAtPosition :: NonEmpty TrainAnchor -> Double -> Maybe TrainAnchor

data LinearExtrapolator

instance Extrapolator LinearExtrapolator where
  -- TODO: this kind of sorting is bullshit; should look up depending on position,
  -- not time (else there's glitches)
  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 $ sortBy (flippedCompare `on` trainAnchorWhen)
    $ NE.filter (\a -> trainAnchorSequence a < positionNow) history
    where difference status = positionNow - (trainAnchorSequence status)
          flippedCompare a b = case compare a b of
            LT -> GT
            GT -> LT
            a  -> a

  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 = unsafePerformIO $ do
  print (scheduledPosition, round $ (expectedProgress - observedProgress) * int2Double (unSeconds expectedTravelTime))
  pure $ (Seconds $ round $ (expectedProgress - observedProgress) * int2Double (unSeconds expectedTravelTime)
         , scheduledPosition)
  where closestPoint =
            minimumBy (compare `on` euclid (trainPingLat, trainPingLong)) line
        nextStop = snd $
            minimumBy (compare `on` fst)
            $ V.filter (\(dist,_) -> dist > 0)
            $ fmap (\stop -> (distanceAlongLine line closestPoint (stationGeopos $ stopStation stop), stop)) tripStops
        lastStop = snd $
            maximumBy (compare `on` fst)
            $ V.filter (\(dist,_) -> dist < 0)
            $ fmap (\stop -> (distanceAlongLine line closestPoint (stationGeopos $ stopStation stop), stop)) tripStops
        line = shapePoints tripShape
        expectedTravelTime =
            toSeconds (stopArrival nextStop) tzseries runningDay
            - toSeconds (stopDeparture lastStop) tzseries runningDay
        expectedProgress = crop $
            seconds2Double (utcToSeconds trainPingTimestamp runningDay
                         - toSeconds (stopDeparture lastStop) tzseries runningDay)
            / seconds2Double expectedTravelTime
            where crop a
                    | a < 0 = 0
                    | a > 1 = 1
                    | otherwise = a
        observedProgress =
            distanceAlongLine line (stationGeopos $ stopStation lastStop) closestPoint
            / distanceAlongLine line (stationGeopos $ stopStation lastStop) (stationGeopos $ stopStation nextStop)
        scheduledPosition =
          (int2Double $ stopSequence lastStop) + observedProgress * (int2Double $ stopSequence nextStop - stopSequence lastStop)

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) = snd
          $ foldl (\(p,a) p' -> (p', a + euclid p p')) (V.head line,0)
          $ V.take (index + 1) line
          where index = fst $ minimumBy (compare `on` (euclid p . snd))
                  $ V.indexed line

-- | convert utc time to seconds on a day, with wrap-around
-- for trains that cross midnight.
utcToSeconds :: UTCTime -> Day -> Seconds
utcToSeconds time day =
  Seconds $ round $ nominalDiffTimeToSeconds $ diffUTCTime time (UTCTime day 0)

secondsNow :: MonadIO m => Day -> m Seconds
secondsNow runningDay = do
  now <- liftIO getCurrentTime
  pure $ utcToSeconds now runningDay

euclid :: Fractional f => (f,f) -> (f,f) -> f
euclid (x1,y1) (x2,y2) = x*x + y*y
  where x = x1 - x2
        y = y1 - y2