1 files changed, 97 insertions, 0 deletions
diff --git a/lib/Extrapolation.hs b/lib/Extrapolation.hs
new file mode 100644
index 0000000..9b3f89f
--- /dev/null
+++ b/lib/Extrapolation.hs
@@ -0,0 +1,97 @@
+{-# LANGUAGE AllowAmbiguousTypes     #-}
+{-# LANGUAGE ConstrainedClassMethods #-}
+{-# LANGUAGE ConstraintKinds         #-}
+{-# LANGUAGE DataKinds               #-}
+{-# LANGUAGE RecordWildCards         #-}
+
+module Extrapolation (Extrapolator(..), LinearExtrapolator, linearDelay) where
+import           Data.Foldable (maximumBy, minimumBy)
+import           Data.Function (on)
+import qualified Data.Map      as M
+import           Data.Time     (Day, UTCTime (UTCTime), diffUTCTime,
+                                nominalDiffTimeToSeconds)
+import qualified Data.Vector   as V
+import           Persist       (Running (..), TrainAnchor (..), TrainPing (..))
+
+import           GHC.Float     (int2Double)
+import           GHC.IO        (unsafePerformIO)
+import           GTFS          (Depth (Deep), GTFS (..), Shape (..), Stop (..),
+                                Time, Trip (..), stationGeopos, toSeconds)
+
+
+
+class Extrapolator a where
+  guessStatusAt :: [TrainAnchor] -> UTCTime -> TrainAnchor
+  guessAnchor :: GTFS -> Running -> TrainPing -> TrainAnchor
+
+data LinearExtrapolator
+
+instance Extrapolator LinearExtrapolator where
+  guessStatusAt history when =
+    minimumBy (compare `on` difference)
+    $ filter (\a -> trainAnchorWhen a > when) history
+    where difference status = diffUTCTime when (trainAnchorWhen status)
+
+  guessAnchor gtfs@GTFS{..} Running{..} ping@TrainPing{..} = TrainAnchor
+    { trainAnchorCreated = trainPingTimestamp
+    , trainAnchorTrip = runningTrip
+    , trainAnchorDay = runningDay
+    , trainAnchorWhen = trainPingTimestamp
+    , trainAnchorDelay = Just (linearDelay gtfs trip ping runningDay)
+    , trainAnchorMsg = Nothing
+    }
+    where Just trip = M.lookup runningTrip trips
+
+linearDelay :: GTFS -> Trip Deep Deep -> TrainPing -> Day -> Int
+linearDelay GTFS{..} trip@Trip{..} TrainPing{..} runningDay = unsafePerformIO $ do
+  print (nextStop, lastStop)
+  print expectedTravelTime
+  -- print (((utcToSeconds trainPingTimestamp runningDay),  toSeconds (stopDeparture lastStop)))
+  print (observedProgress, expectedProgress)
+  pure $ round $ (expectedProgress - observedProgress) * int2Double expectedTravelTime
+  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 trainPingTimestamp
+            - toSeconds (stopDeparture lastStop) tzseries trainPingTimestamp
+        expectedProgress =
+            (int2Double ((utcToSeconds trainPingTimestamp runningDay)
+                         - toSeconds (stopDeparture lastStop) tzseries trainPingTimestamp))
+            / int2Double 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)
+
+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 -> Int
+utcToSeconds time day =
+  round $ nominalDiffTimeToSeconds $ diffUTCTime time (UTCTime day 0)
+
+euclid :: Fractional f => (f,f) -> (f,f) -> f
+euclid (x1,y1) (x2,y2) = x*x + y*y
+  where x = x1 - x2
+        y = y1 - y2