`frequencies_to_stop_times` slow execution speed in large datasets

[cseveren]

Thanks for the awesome and helpful package! I get a lot of value out of frequencies_to_stop_times, but I've noticed that it's quite slow when working at scale (i.e., it takes several hours on my reasonably powerful desktop to go from 37k stop times and 1k trips to 4.4mil stop times and 140k trips for recent Mexico City freq-based GTFS feeds). This is almost assuredly due to the use of rbind here to append row by row onto the eventually quite huge stop_times data.table which is used to store and eventually output results:

 stop_times <- rbind (stop_times, stop_times_trip_i)

I propose three alternatives, which I've tried to implement locally on my machine but I've not worked with rcpp before and so I couldn't make work. I can't quite tell if stop_times_trip_i is always 1 row long or not. If so:

• Preallocate stop_times with some upper bound number of rows (perhaps user set), and then simply replace the above with:

   stop_times[i, ] <- stop_times_trip_i
  

Else (if it's number of rows vary), here are two other ideas:

• Preallocate stop_times with some upper bound number of rows (perhaps user set), and replace a group of rows, as above
• Accumulate results in smaller data.tables and then rbind, there's an opportunity for a couple of levels of recursion. Specifically, the for trips block could be (sorry for the stupid names):

    for (trip in trips)  {

        # n is to  be added to trip_id
        n <- 1

        # order frequencies table by trip_id and start_time
        frequencies_trip <-
            gtfs_cp$frequencies [order (trip_id, start_time)] [trip_id == trip]

        # in case end_time of the previous period and start_time of the next are
        # equal:
        frequencies_trip [end_time ==
                          data.table::shift (start_time, 1,
                                             type = "lead"),
                            end_time := end_time - 1]

        stop_times_trip <-
            gtfs_cp$stop_times [order (stop_sequence)] [trip_id == trip]

        # in case of the first arrival time > 0, then reset them:
        if (stop_times_trip [1] [["arrival_time"]] > 0)  {
            stop_times_trip [, c ("arrival_time", "departure_time") :=
                             list(
                arrival_time - stop_times_trip [1] [["arrival_time"]],
                departure_time - stop_times_trip [1] [["arrival_time"]]
                )]
        }

        start_t <- min (frequencies_trip$start_time)

        headway <- headway_old <- frequencies_trip [1] [["headway_secs"]]

        ## NEW LINE ##
        stop_times_small <- stop_times

        for (i in row (frequencies_trip))  {

            end_t <- frequencies_trip [i] [["end_time"]]
            headway <- frequencies_trip [i] [["headway_secs"]]

            # in order to ensure a 'smooth' transition between frequency periods
            ifelse (
                headway_old - start_t +
                    frequencies_trip [i] [["start_time"]] < headway,
                start_t <- start_t - headway_old + headway,
                start_t <- frequencies_trip [i] [["start_time"]]
            )

            ## NEW LINE ##
            stop_times_small_small <- stop_times

            # multiply stop_times for all trips based on a given frequency
            while (start_t < end_t)  {

                stop_times_trip_i <-
                    data.table::copy (stop_times_trip) [, c ("arrival_time",
                                                             "departure_time",
                                                             "trip_id_f") :=
                            list ((arrival_time + start_t),
                                  (departure_time + start_t),
                                  paste (trip_id, n, sep = "_"))]

                n <- n + 1

                # EDITED BELOW #
                stop_times_small_small <- rbind (stop_times_small_small, stop_times_trip_i)

                start_t <- start_t + headway
            }
            # EDITED BELOW #
            stop_times_small <- rbind (stop_times_small, stop_times_small_small)

            headway_old <- headway

        }
        # EDITED BELOW #
        stop_times <- rbind (stop_times, stop_times_small_small)

    }

This last option, while not totally efficient, will probably substantially increase performance without needing to explicitly deal with preallocation. Again, thanks for the great package!

[cseveren]

Ok -- In my use case, following the third option above (not exactly, I made a couple errors) sped execution time ~100x. It still drags toward the end bc it's accumulating to a progressively larger and larger data.table, but it's enough of a boost to make my use case work.

My fork is here: https://github.com/cseveren/gtfs-router - Let me know if you'd like me to submit a pull request. NB: I also added display tracker/counter on original trips.

[dhersz]

Bringing the discussion that arose in r5r's repo to this issue:

The development version of gtfstools also includes a frequencies_to_stop_times() function. It's currently a bit faster than gtfsrouter's:

path <- system.file("extdata/spo_gtfs.zip", package = "gtfstools")

gtfstools_gtfs <- gtfstools::read_gtfs(path)
gtfsrouter_gtfs <- gtfsrouter::extract_gtfs(path, quiet = TRUE)

microbenchmark::microbenchmark(
  cvt_gtfstools_gtfs <- gtfstools::frequencies_to_stop_times(gtfstools_gtfs),
  cvt_gtfsrouter_gtfs <- gtfsrouter::frequencies_to_stop_times(gtfsrouter_gtfs),
  times = 5L
)
#> Unit: milliseconds
#>                                                                           expr
#>     cvt_gtfstools_gtfs <- gtfstools::frequencies_to_stop_times(gtfstools_gtfs)
#>  cvt_gtfsrouter_gtfs <- gtfsrouter::frequencies_to_stop_times(gtfsrouter_gtfs)
#>         min        lq      mean   median         uq        max neval
#>    777.2107   793.785   875.153   925.58   937.2485   941.9406     5
#>  26631.9393 31469.455 31593.529 32363.29 32884.5130 34618.4487     5

Also, there seems to be a bug in gtfsrouter's, since it doesn't update the frequencies table:

cvt_gtfstools_gtfs$frequencies
#> NULL

cvt_gtfsrouter_gtfs$frequencies
#>         trip_id start_time end_time headway_secs
#>   1: CPTM L07-0      14400    17940          720
#>   2: CPTM L07-0      18000    21540          360
#>   3: CPTM L07-0      21600    25140          360
#>   4: CPTM L07-0      25200    28740          360
#>   5: CPTM L07-0      28800    32340          360
#>  ---                                            
#> 700:  5290-10-1      79200    82740         1200
#> 701:  5290-10-1      82800    86340         1200
#> 702:  6450-51-0      18000    21540         3600
#> 703:  6450-51-0      21600    25140         3600
#> 704:  6450-51-0      25200    28740         3600

It would be interesting to hear from @cseveren if gtfstools' function works fine for his use case. I could also work on a PR if @mpadge is happy with my implementation and wants to use it.

P.S.: I don't want to sound arrogant/snob coming to this issue and showing that "my function" is faster, or anything like that. Before implementing frequencies_to_stop_times() in gtfstools I admittedly looked at gtfsrouter implementation and took a lot of inspiration from it. This conversation randomly came up in r5r as I was looking at some old open issues, and I decided that moving it to this issue would be beneficial because I'd be able to properly report what seems to me like a bug and we would be able to discuss both implementations in further details.

[cseveren]

Much faster than my hack-y solution! Took only 20-30 seconds.

[mpadge]

Wow, thanks @cseveren for raising this important issue, and especially thanks to @dhersz for chiming in and offering a clearly superior solution, for which you definitely do not sound arrogant 😆 I was originally going to just ask for a PR of the gtfstools function into here (mostly because I want to keep this package as light on dependencies as possible, and importing gtfstools would then brings in the ultimately heavy dep of sf, and a few others). But then I had a bit of a look and realised it might help all of us if I try to rewrite the function in C++. The version here was just assembled as a hack, and the slowness is precisely for reasons you asserted @cseveren, but that's something that could be totally circumvented by porting the function to C++.

My suggestion would then be for me to start doing that in a branch here which we can use to compare with the current gtfstools implementation. I'll make sure I do it cpp11-style, rather than Rcpp, so it can easily be plugged straight in there if it turns out to offer significant speed gains. Only catch: I'm about to be away from normal routine for a few weeks, so may have little opportunity before March (or later) to get cracking on this. Hope that's okay!

[mpadge]

Those commits finally finish converting the "frequencies_to_stop_times" routine to C++. It's slightly faster than @dhersz's gtfstools version, although not by a huge amount (around 20%). I think the big advantage is that this operation is much clearer when coded in C++, and so actually takes much less code as well - around 70 lines of very clear C++ plus 78 lines of R, most of which are input checks. This compares with gtfstools R function of 336 lines, and former function here of 170 lines that was (i) confusing, and (ii) didn't work properly.

@dhersz Our two packages do not give identical results, but largely because this procedure requires a few seemingly arbitrary decisions to be made. These notably include:

1. What to to if "stop_time" for a "frequencies" entry is marginally less than some integral of "start_time + n * headway" - this does happen, and I've opted to include an addiitonal trip in those cases.
2. What do to when "stop_times" entries for "trip_id" values which are in the frequencies table have "arrival/departure_time" values of "0", yet "trip_id" values flagged in "frequencies" tables with "exact_times" values of 1. These latter values of 1 are supposed to indicate exact arrival/departure times, which should be used as given, but obviously arrival/departure time values of 0 are nonsense. I've opted to ignore "exact_times" flags in "frequencies" tables, becuase these seem to be frequently erroneous.

[mpadge]

Re-opening because there is one pronounced negative effect on other routines. The routine implemented here appends _[0-9]+ values to the end of trip_id entries in the "stop_times" table. The internal routine used to construct timetables relies on matching "trip_id" values, and was modified here to grep matches via an lapply call, because there are commonly too many trip_id values to form a single compound grep pattern. This lapply(ptn, grep) call is, however, really slow, so needs to be replaced with a better way of mapping potentially modified trip_id values back to original, unmodified forms.