LASPATED: A Library for the Analysis of Spatio-Temporal Discrete Data

2.1 Model without Covariates

Each time interval $t\in \mathcal{T}$ has a duration ${\mathcal{D}}_t$ (in time units). It is assumed that data are observed for multiple occurrences of the same interval t, and that each time data are observed for a type $c\in \mathcal{C}$, zone $i\in \mathcal{I}$, and interval $t\in \mathcal{T}$, all such space-time points in the same time interval are recorded in the data. For example, suppose that $\mathcal{T}$ forms a partition of the week in time intervals of 1 hour each. Then it is assumed that the data contain observations of emergencies of type c in zone i during hour t for multiple weeks, and that each time emergency arrivals of type c in zone i are observed during a particular hour t of a particular week, all the emergency arrivals of the same type in the same zone during the same hour t of that week are recorded. (It is planned to make provision for censored data in future work.) All the emergency arrivals of a particular type c in a particular zone i during a particular hour t of a week together are called an observation. For each $c\in \mathcal{C}$, $i\in \mathcal{I}$, and $t\in \mathcal{T}$, let N_c,i,t denote the number of observations for type c, zone i, and time interval t, and let these observations be indexed by $n\in {\mathcal{N}}_{\mathit{\text{c,i,t}}}:=\{1,\mathrm{\dots },N_{\mathit{\text{c,i,t}}}\}$. For each $c\in \mathcal{C}$, $i\in \mathcal{I}$, $t\in \mathcal{T}$, and $n\in {\mathcal{N}}_{\mathit{\text{c,i,t}}}$, let M_c,i,t,n denote the number of points (arrivals) for observation n of type c, zone i, and time interval t, and let $M_{\mathit{\text{c,i,t}}}:={\sum }_{n=1}^{N_{\mathit{\text{c,i,t}}}}{M}_{c,i,t,n}$ denote the total number of points over all observations for type c, zone i, and time interval t.

Assume that $\{M_{c,i,t,n}:c\in \mathcal{C},i\in \mathcal{I},t\in \mathcal{T},n\in {\mathcal{N}}_{\mathit{\text{c,i,t}}}\}$ are independent (but not necessarily identical) Poisson distributed random variables. Let ${\lambda }_{\mathit{\text{c,i,t}}}$ denote the mean number of points per length of time (such as per hour) for type c, zone i, and time interval t. Then random variable M_c,i,t,n is Poisson distributed with mean ${\lambda }_{\mathit{\text{c,i,t}}}{\mathcal{D}}_t$. Let $\lambda :=({\lambda }_{\mathit{\text{c,i,t}}},c\in \mathcal{C},i\in \mathcal{I},t\in \mathcal{T})$. Then the likelihood function is

and intensities $\lambda$ that maximize the log-likelihood are the same intensities that solve

A typical issue with such applications is that the distribution of data is far from uniform — there are a few combinations of type c, zone i, and time interval t with many observations, but for most combinations of type c, zone i, and time interval t there are very few observations. In such cases, it may be advantageous to use data from “neighboring” observations, or to bring additional data to bear on the estimation problem.

For example, suppose that $\mathcal{T}$ forms a partition of the week in time intervals of 1 hour each, but that it is expected that many hours of the week are similar to other hours of the week in terms of arrival rates. A simple approach to incorporate such an idea is to partition $\mathcal{T}$ into a collection $\mathcal{G}$ of subsets of $\mathcal{T}$ in such a way that it may be reasonable to expect that, for each $c\in \mathcal{C}$, $i\in \mathcal{I}$, and $G\in \mathcal{G}$, the values of ${\lambda }_{\mathit{\text{c,i,t}}}$ for different $t\in G$ will be close to each other (but not necessarily the same). Let $W_G\ge 0$ denote a similarity weight for $G\in \mathcal{G}$. Another approach to incorporate such an idea is to specify a similarity weight $w_{t,t^{’}}\ge 0$ for each (unordered) pair $t,t^{’}\in \mathcal{T}$. For example, $w_{t,t^{’}}=w>0$ if t and t’ are neighboring time intervals, and $w_{t,t^{’}}=0$ otherwise. Similarly, it may be expected that many zones are similar to other zones in terms of arrival rates. Such an idea can also be incorporated by partitioning the set $\mathcal{I}$ of zones into subsets, each with its own similarity weight, or by specifying for each pair $\mathit{\text{i,t}}\in \mathcal{I}$, a similarity weight $w_{\mathit{\text{i,t}}}\ge 0$. An example loss function with similarity regularization that uses the first approach for time intervals and the second approach for zones is given by

Example 2.1: Consider a simple example in which $\mathcal{I}=\{0,\mathrm{\dots },I\}$,  and $\mathcal{C}$ and $\mathcal{T}$ are singletons, that is, $\lambda$ is estimated as a function of space only. Suppose that the number of observations $N_i=N$ is the same for all $i\in \mathcal{I}$. Note that the empirical estimates, without smoothing, are ${\widehat{\lambda }}_i^E=M_i/N$. For smoothing, choose $w_{0,i}=1$ for all $i\in \mathcal{I}$, and $w_{\mathit{\text{i,t}}}=0$ otherwise, that is, zone 0 is a neighbor of all other zones, but none of the other zones are neighbors of each other. Then the loss function $\mathrm{\ell }$ in (6) reduces to

It follows from the first-order optimality conditions that

Unless ${\lambda }_i={\lambda }_0$ for all i, it follows from $N{\sum }_{i\in \mathcal{I}\setminus \{0\}}{\left({\lambda }_i–{\lambda }_0\right)}^2>0$ that ${\sum }_{i\in \mathcal{I}}{\lambda }_i<{\sum }_{i\in \mathcal{I}}{\widehat{\lambda }}_i^E$, that is, the estimated total rate is biased low.

To counter the bias illustrated in Example 2.1, one may want to use the data to constrain the estimates $\lambda$. For example, taking duration ${\mathcal{D}}_t$ of time interval t into account, one may want to constrain the estimates $\lambda$ so that for every type $c\in \mathcal{C}$, the sum of the intensities ${\lambda }_{\mathit{\text{c,i,t}}}{\mathcal{D}}_t$ is equal to the sum of the empirical intensities ${\widehat{\lambda }}_{\mathit{\text{c,i,t}}}^E{\mathcal{D}}_t=M_{\mathit{\text{c,i,t}}}/N_{\mathit{\text{c,i,t}}}$ over zones $i\in \mathcal{I}$ and time intervals $t\in \mathcal{T}$. The corresponding constraint is

Intensity estimates $\lambda$ are then obtained by solving

2.2 Model with Covariates

Each type c, zone i, and time interval t may have covariates that are correlated with the arrival rates ${\lambda }_{\mathit{\text{c,i,t}}}$, and data of these covariates can be used to partly compensate for sparse data. For example, emergency arrival rates in different zones and time intervals can be expected to be correlated with the population and other measures of economic activity in the zones, as well as with festivals and other events during the time intervals. For each $c\in \mathcal{C}$, $i\in \mathcal{I}$, and $t\in \mathcal{T}$, let $x_{\mathit{\text{c,i,t}}}:=(x_{c,i,t,1},\mathrm{\dots },x_{c,i,t,K})$ denote the covariate values of type c, zone i, and time interval t. For example, $x_{c,i,t,1}$ may be the population count with home addresses in a zone i, and $x_{c,i,t,2}$ may be an indicator that a major sports event is scheduled in a zone i during a time interval t. Then consider the model

where $\beta =({\beta }_1,\mathrm{\dots },{\beta }_K)$ are the model parameters. Let ${\mathcal{X}}_{\mathit{\text{c,i,t}}}$ denote the set of all possible values of x_c,i,t. Often ${\mathcal{X}}_{\mathit{\text{c,i,t}}}$ can be chosen to be a polyhedron. Note that it should hold that

To facilitate such a model, let $N={\sum }_{c\in \mathcal{C}}{\sum }_{i\in \mathcal{I}}{\sum }_{t\in \mathcal{T}}{N}_{\mathit{\text{c,i,t}}}$ denote the total number of observations, and let these observations be indexed $n=1,\mathrm{\dots },N$. For each observation $n\in \{1,\mathrm{\dots },N\}$, let M_n denote the number of arrival points for observation n and let $x^n:=(x_1^n,\mathrm{\dots },x_K^n)$ denote the covariate values of observation n. Then the negative log-likelihood function is given by

As for the model without covariates, one may want to use the data to constrain the estimates $\lambda$. For example, one may add the constraint that for every type $c\in \mathcal{C}$, the sum over all zones and all time intervals of the estimated rates is equal to the sum over all zones and all time intervals of the empirical rates, that is the constraint

Next, we provide two examples of such models.

Example 2.2: Index $\mathcal{C}=\{1,\mathrm{\dots },|\mathcal{C}|\}$, $\mathcal{T}=\{1,\mathrm{\dots },|\mathcal{T}|\}$, and let $K_1:=|\mathcal{C}\times \mathcal{T}|$. For each $c\in \mathcal{C}$, $t\in \mathcal{T}$, and $k=(c–1)|\mathcal{T}|+t$, let $x_k^n$ be the number of people resident (population count) in the zone of observation n if observation n is for type c and time interval t, and $x_k^n=0$ otherwise. The next covariate is a set of occupational land use areas (in km²), for instance the areas of commercial activities and public facilities, of industrial activities, and of undeveloped land. Index the occupational land uses by $\mathcal{O}=\{1,\mathrm{\dots },|\mathcal{O}|\}$. For each $c\in \mathcal{C}$, $t\in \mathcal{T}$, $m\in \mathcal{O}$, and $k=K_1+(c–1)|\mathcal{T}||\mathcal{O}|+(t–1)|\mathcal{O}|+m$, let $x_k^n$ be the area of occupational land use m in the zone of observation n if observation n is for type c and time interval t, and $x_k^n=0$ otherwise. Then ${\sum }_{k=1}^{K_1} {\beta }_k{x}_k^n$ denotes the forecasted number of arrivals for the type of observation n in the zone of observation n during the time interval of observation n due to people being in the residential area, which is modeled as proportional to the number of people resident in the zone of observation n with a proportionality coefficient that depends on the type and the time interval. Similarly, ${\sum }_{k=K_1+1}^{K_1(1+|\mathcal{O}|)}{\beta }_k{x}_k^n$ denotes the forecasted number of arrivals for the type of observation n in the zone of observation n during the time interval of observation n due to people being in the different occupational areas which is modeled as proportional to the areas of occupational land use in the zone of observation n with a proportionality coefficient that depends on the type and the time interval. If $K_1|\mathcal{O}|$ is large, the number of parameters ${\beta }_k$ of the model specified in this example is also large. Similar to (6), similarity regularization can be used to estimate the parameters, for example, by minimizing a loss function such as

Example 2.3: This example will be used to demonstrate the modeling of arrivals of emergency calls to an emergency medical service. The sets $\mathcal{C}$, $\mathcal{I}$, and $\mathcal{O}$ are the same as in Example 2.2. Let $\mathcal{T}$ denote the indices of discrete time periods during a day, for example, $\mathcal{T}=\{1,\mathrm{\dots },48\}$ if each day is discretized into 30-minute intervals. Let $\mathcal{D}$ denote the set of indices of the (normal) days of the week, as well as indices for special days such as holidays. Thus, the cardinality of $\mathcal{D}$ is 7 plus the number of special days. A pair $(d,t)\in \mathcal{D}\times \mathcal{T}$ specifies a time interval. For each zone $i\in \mathcal{I}$, consider covariates $x_i\in {\mathbb{R}}^{1+|\mathcal{O}|}$, where $x_i(1)$ is the population count in zone i, and $x_i(j+1)$ is the area (in km²) of occupational land use j in zone i for $j=1,\mathrm{\dots },|\mathcal{O}|$. Here, we assume that these data (population and land type areas) do not depend on time, which is reasonable for moderate time periods. Then, for each type c, zone i, day d, and time period t, the arrival intensity is given by ${\lambda }_{c,i,d,t}{\mathcal{D}}_t={\beta }_{c,d,t}^{\top }{x}_i$ for some vector ${\beta }_{c,d,t}\in {\mathbb{R}}^{1+|\mathcal{O}|}$). Furthermore, for each type c, zone i, day d, and time period t, let N_c,i,d,t denote the number of observations, let M_c,i,d,t,n denote the number of arrival points for observation n, and let $M_{c,i,d,t}:={\sum }_{n=1}^{N_{c,i,d,t}}{M}_{c,i,d,t,n}$. Then the negative of the log-likelihood function is given by

Note that call rates ${\beta }_{c,d,t}^{\top }{x}_i$ should be positive. Additionally, for each type c, day d, and time period t, the coefficients ${\beta }_{c,d,t}(1)$ that represent the ratio of emergencies to population should be small (say less than 1), which implies the constraints $0\le {\beta }_{c,d,t}(1)\le 1$. Therefore, the estimation problem is to solve the optimization problem

Depending on the solver, better numerical performance may be obtained by replacing the constraints of (14) with

for some $\epsilon >0$ sufficiently small.

2.3 Solving the Optimization Problems

The calibration of models such as those of Section 2.1 requires solving an optimization problem (8) of the form $\min \{\mathrm{\ell }(\lambda ):\lambda >0\}$. The problem can also be formulated as

where C is a closed convex set such as

for a chosen $\epsilon >0$ and e a vector of ones. The objective function $\mathrm{\ell }$ of (16) is differentiable on such C.

The calibration of models such as those of Section 2.2 requires solving an optimization problem (14) of the form

where B is given by constraints such as (15). In both cases, the problem is a convex optimization problem with a differentiable objective function on the feasible set.

The algorithms described below for solving problems such as (16) take regularization weights such as W_G and w_i,j as input. LASPATED also provides cross validation functions, that facilitate use of the data to select regularization weights.

Projected gradient with line search. LASPATED provides the following two methods for solving problems (16) and (18): (1) a projected gradient method with line search along a feasible direction, and (2) a projected gradient method with line search along the boundary. These methods follow (25) except for the update of $\mathrm{\Delta }$. Here, we provide the pseudocode for these methods for a problem of form (16). Let ${\mathrm{\Pi }}_C(\lambda )$ denote the projection of the point $\lambda$ onto a closed convex set C. The pseudocode presents methods that are run for a chosen number of iterations and that start from any feasible ${\lambda }_0\in C$. For alternative stopping criteria, see the discussion at the end of this section.

Note that for problem (16) with $C=\{\lambda :\lambda \ge \epsilon \mathbf{\text{e}}\}$, the projection ${\mathrm{\Pi }}_C(\lambda )=\max \{\lambda ,\epsilon \mathbf{\text{e}}\}$ onto C is easy to compute, where the max is taken componentwise. Furthermore, the first derivatives of $\mathrm{\ell }$ are given by

where G(t) denotes the group of time intervals to which t belongs. For problem (18), the first derivatives of $\mathcal{L}$ are given by

For problem (18) with $B=\{\beta :{\beta }_{c,d,t}^{\top }{x}_i\ge \epsilon ,0\le {\beta }_{c,d,t}(1)\le 1,\forall c\in \mathcal{C},d\in \mathcal{D},t\in \mathcal{T}\}$, the projection ${\mathrm{\Pi }}_B(\beta )$ onto B is computed by solving a convex quadratic problem.

Stopping test for projected gradient method with line search. LASPATED makes provision for several stopping criteria such as convergence of the objective function values, a maximum number of iterations, or stopping when optimality conditions are approximately satisfied. We give more details about this latter criterion. Consider an optimization problem

where $\mathrm{\ell }$ and C are convex, and C is compact. At the end of iteration k, points ${\lambda }_0,\mathrm{\dots },{\lambda }_k$ have been generated, and the upper bound

on ${\mathrm{\ell }}_{\ast }$ can be computed. One can minimize the linear approximation of $\mathrm{\ell }$ at ${\lambda }_k$ to compute the lower bound

or one can minimize the piecewise linear approximation of $\mathrm{\ell }$ at ${\lambda }_0,\mathrm{\dots },{\lambda }_k$ to compute the lower bound

The algorithm is stopped when

thus obtaining an $\epsilon$-optimal solution.

3.1 Periodic with Equal Length Time Intervals

The simplest time discretization in LASPATED uses a periodic pattern, in which the duration of the periodic pattern is partitioned into time intervals of equal length. That is, each subset of the time discretization is a single time interval, and all these time intervals have the same length.

3.2 Periodic with Unequal Length Time Intervals

Another time discretization in LASPATED also uses a periodic pattern, in which the duration of the periodic pattern is also partitioned into time intervals, but the time intervals may have unequal lengths. As in the previous method, each subset of the time discretization is a single time interval, but unlike the previous method, all these time intervals do not have the same length.

3.3 Customized Subsets

LASPATED makes provision for customized subsets to facilitate holidays and special events. Each customized subset is assigned a unique index number, and consists of one or more time intervals. Each time interval is specified by its start time and its end time, as well as the index of the subset that the time interval belongs to. Time points that do not belong to any customized interval, belong to the customized subset with index 0. LASPATED also allows time intervals to repeat. For example, a time interval that starts and ends on date 2016-01-01 (New Year’s Day) may be specified to repeat each year. If the time interval is specified as repeating yearly, then for all observations on the same day of the year (such as January 1), the time intervals’ start time during the day, end time during the day, as well as the index of the subset that the time intervals belong to, will be the same for all years.

4.1 Defining Borders

The first step for space discretization is to choose a coordinate system and to specify the border of a region that contains the locations of all the points in the dataset. The border specifies the region that will be discretized in space. The border can be specified with LASPATED using different methods.

4.2 Equal Sized Rectangular Space Discretization

The simplest space discretization in LASPATED partitions the region into equal sized rectangles, in such a way that adjacent rectangles share a common face, that is, if the boundaries of rectangle A and rectangle B intersect, then either rectangle A and rectangle B intersect in one (corner) vertex, or rectangle A and rectangle B share an edge between 2 vertices. For example, the rectangles $A=[0,2]\times [0,2]$ and $B=[2,4]\times [1,3]$ are adjacent but do not share a common face. If the region is not a union of these equal sized rectangles, then the intersections of the rectangles with the region may not be equal sized. For example, Figure 1 shows a discretization of a custom region containing the city of Rio de Janeiro into $10\times 10=100$ rectangles, and Figure 7 shows a discretization of the same region into $100\times 100=10000$ rectangles. Both discretizations were obtained with LASPATED.

Figure 7

4.3 Equal Sized Hexagonal Space Discretization

Another simple space discretization in LASPATED partitions the region into equal sized hexagons, in such a way that adjacent hexagons share a common edge. (The Uber Python package H3 was used for the discretization.) Figure 2 shows a discretization of a region containing the city of Rio de Janeiro into hexagons using a scale parameter equal to 7, and Figure 8 shows a discretization of the same region into hexagons using a scale parameter equal to 8. Both discretizations were obtained with LASPATED.

Figure 8

4.4 Customized Space Discretization

LASPATED also facilitates space discretization with customized subregions. For example, Figure 3 obtained with LASPATED displays a customized discretization of the city of Rio de Janeiro into 160 administrative districts.

4.5 Discretization using Voronoi diagrams

LASPATED also provides space discretization with Voronoi diagrams. Figure 9 obtained with LASPATED displays a discretization based on the Voronoi diagram given by the locations of ambulance stations in Rio de Janeiro. Each subregion includes the set of points that are closest to a specific station.

Figure 9

6.1 Compilation and Dependencies

C++ calibration functions use the following libraries:

• Gurobi (https://www.gurobi.com), an optimization tooblox used only in the model with covariates;

• xtensor and xtl libraries.

When compiling the code, the user needs to pass the location of the header files for LASPATED, xtl, xtensor and gurobi, as well as the location for Gurobi libraries and their respective flags. xtl and xtensor require the C++14 standard and its flag must also be passed. The example below compiles an executable called laspated from a main.cpp that includes laspated.h, and uses the Gurobi solver.

$ g++ -o laspated main.cpp -std=c++14
    -I<gurobi_header_dir>
    -I<laspated_header_dir>
    -I<xtl_header_dir>
    -I<xtensor_header_dir>
    -L<gurobi_libs_dir>-lgurobi_c++ -lgurobi100

If Gurobi is not used, the calibration functions for the model without covariates are used by setting the flag -DUSE_GUROBI equal to zero. The following command compiles the code in such a case:

$ g++ -o laspated main.cpp -std=c++14
    -DUSE_GUROBI=0 -I<laspated_header_dir>
    -I<xtl_header_dir> -I<xtensor_header_dir>

We provide five test models in the source files test_problems.cpp and test_problems.h in the Model_Calibration directory. They were used to assert the quality of the projected gradient algorithm.

6.2 LASPATED C++ Application

We also provide a LASPATED binary app that allows running the calibration functions. The user must provide a set of files describing the model and a set of configuration options, either via a configuration file or command line. The options describe the algorithmic parameters, model attributes, and file locations. The valid options are the following:

• double EPS: tolerance used in the projection algorithm;

• double sigma: parameter used in the line search of the projected gradient method;

• double accuracy: gap used as stopping criterion for the projected gradient method;

• int max_iter: maximum number of iterations for the projected gradient method;

• double lower_lambda: lower bound for decision variables in the projected gradient method;

• double upper_lambda: upper bound for decision variables in the projected gradient method;

• double beta_bar: initial step size for the projected gradient method;

• double cv_proportion: proportion of samples used in each iteration of cross validation;

• output_file: path where the LASPATED app will save the estimated intensities.

• string model_type: the type of model being solved, must be either “reg” (for the model with covariates) or “no_reg” (for the model without covariates);

• string method: the type of method being used. For the model with covariates, only “calibration” is supported (running the projected gradient algorithm once). For the model without covariates, can be either “calibration” or “cross_validation”;

• string algorithm: the algorithm used. Can be either “feasible” for Armijo search along the feasible direction or “boundary” for Armijo search along the boundary. However “cross_validation” only supports the “feasible” option;

• string info_file: path of file containing general information about the model. See below for a precise description of the corresponding file.

• arrivals_file: path of file containing the samples of arrivals. See below for a precise description of the corresponding file.

• neighbors_file: path of file containing the description of each zone, including its type, its covariates, and its neighbors. See below for a precise description of the corresponding file.

• durations_file: path of file containing the durations for each time index. See below for a precise description of the corresponding file.

• alpha_regions_file: path of file containing the space penalizations $R\times R$ matrix. See below for a precise description of the corresponding file.

• time_groups_file: path of a file containing a description of time groups and the corresponding penalties in the regularized model from Section 2.1. See below for a precise description of the corresponding file.

• double duration: array with the duration of each period.

• cv_weights_file: path of a file containing the weights to be used in cross validation. See below for a precise description of the corresponding file.

• relax_empirical_fix: flag that, when set, relaxes the constraint that makes the sum of estimated rates per class equal to the empirical rate by class.

Note that to use the model with covariates, the user only needs to specify the info, arrivals, neighbors, and duration files. Additionally, for the model without regressors, the user must specify the files alpha_regions_file and time_groups_file. When performing cross validation, the file cv_weights_file is also needed.

The file info_file must contain two lines. The first line must contain, in that order, separated by spaces, the number of periods $T$ per day, the number of week days $G$, the number of zones $R$, the number of arrival types $C$, the number of covariates $J$ and the number of holidays $H$ for each year. The second line contains, for each day of the week and holiday, the number of observations $N_g$ for that day or holiday. Values of index $g$ from 0 to 6 are used for days of the week while holidays have indexes in $\{7,8,\dots ,7+H\}$.

The file arrivals_file contains the samples observed, where each line describes a sample. Each line contains, in that order, separated by spaces, the period $t\in [0,T)$, the day $g\in [0,G)$, the zone $r\in [0,R)$, the arrival type $c\in [0,C)$, the sample index $j\in [0,N_g)$, the number of arrivals observed, and a flag indicating whether the day is a holiday or not.

The file neighbors_file contains the description of each zone. Each line contains, in that order, the zone index, two values representing the coordinates of the zone centroid, a number indicating the zone’s land type, J values corresponding to the values of all J covariates, followed by a set of pairs of values (index, distance) indicating the distances of the considered zone to its neighbors (identified by their indexes).

The file alpha_regions_file describes the space penalization parameters. It contains a $R\times R$ matrix (with entries separated by spaces) where entry $(i,j)$ corresponds to the weight in the objective function of the model without covariates if zones $i$ and $j$ are neighbors.

The file time_groups_file contains the description of time groups and their corresponding penalties. The first line must have the number $|\mathcal{G}|$ of time groups. The subsequent $D\ast T$ lines must have the group index of each period. Then the subsequent $\mathcal{G}$ lines may contain the weight of each group. If the method is set to “cross_validation”, the reading of time weights is skipped.

File cv_weights_file contains the weights to be used in cross validation, separated by spaces. Each weight must be non-negative.

The durations file contains the duration of each time index, in hours. The first line contains, separated by spaces, the number of day indexes and the number of time indexes.

The LASPATED app is only supported on Linux operating systems. The user must have the GCC compiler and the Boost library. In order to use the model with covariates, the user must have Gurobi installed. The directory Model_Calibration/Cpp contains a Makefile that is used to compile the app with:

$ make GUROBI_LIB_VER= -lgurobi110
    GUROBI_HOME=/opt/gurobi1101/linux64
INCLUDE_BOOST=/usr/include/boost

The Makefile uses GUROBI_HOME to find Gurobi directory, GUROBI_LIB_VER to find the Gurobi C library version and INCLUDE_BOOST to find the Boost directory. The values displayed above are actually the default values set in the Makefile, so if Boost and Gurobi are installed in those locations, then the make command will compile the code. Once LASPATED app is successfully compiled, it can be run with command:

$ ./laspated --model_type=no_reg
    --method=calibration
    --info_file=info.txt
    --arrivals_file=arrivals.txt
    --neighbors_file=neighbors.txt
    --alpha_regions_file=alpha_regions.txt
    --time_groups_file=time_groups.txt
    --durations_file=durations.txt
    --output_file=output.txt

to calibrate the model without covariates with space penalizations given in file alpha_regions.txt (see above) and penalization between time groups given in file time_groups.txt (see above). The command

$ ./laspated
    --model_type=no_reg
    --method=cross_validation
    --info_file=info.txt
    --arrivals_file=arrivals.txt
    --neighbors_file=neighbors.txt
    --time_groups_file=time_groups.txt
    --durations_file=durations.txt
    --cv_weights_file=cv_weights.txt
    --output_file=output.txt

uses cross validation to calibrate the model without covariates (and therefore computes penalizations). Finally, the command

$ ./laspated
    --model_type=reg
    --method=calibration
    --info_file=info.txt
    --arrivals_file=arrivals.txt
    --neighbors_file=neighbors.txt
    --durations_file=durations.txt
    --output_file=output.txt

calibrates the model with covariates.

More precisely, the first example performs the model estimation for the model without covariates (option model_type=no_reg) and saves the model parameters in output.txt (option output_file). The second example performs the model estimation without covariates using cross validation (option method=cross_validation). The third example calibrates the model with covariates (model_type=reg). Users can also provide the options via configuration file with the option -f.

After running the LASPATED app, the output is saved in the path provided by option output_file. If the model without covariates is used, the output contains in each line the indexes $c\in [0,C)$, $r\in [0,R)$, $t\in [0,D\ast T)$, and the intensity ${\lambda }_{c,r,t}$ estimated by LASPATED. If the model with covariates is used, then each line correspond to indexes $c\in [0,C)$, $d\in [0,D)$, $t\in [0,T)$, $j\in [0,J)$, and the parameter ${\beta }_{c,d,t,j}$ estimated by LASPATED.

6.3 Demo Scripts

LASPATED repository contains a directory with demo scripts for the main functionalities that can be run with

python demo_discretization.py
    <border> <geo_disc>

to test discretization functions.

The option <border> must be one of the following border types: rectangle, convex and custom. The option <geo_disc> must be one of the discretization types: rectangle, hexagon and custom. This demo creates, based on our medical emergencies’ dataset, the discretization <geo_disc> of the region defined by the border <border>, and saves the files that can be used by the model estimation application. After compiling the C++ application, the following demos run the model estimation functions:

python demo_calibration_regularized.py
python demo_calibration_covariates.py

These scripts create the model files and invoke the C++ application that performs the model estimation. The model parameters are saved into files in the calib_data_covariates and calib_data_regularized directories.

6.4 Replication Script

We also provide, in the repository https://github.com/vguigues/Laspated_Replication, a python script that uses LASPATED to reproduce the results from the Quality Control and Reuse Potential sections. Before running the script, the first step is to compile the C++ code. If you have a GCC compiler installed, this can be done by running the following commands from the Laspated_Replication root directory:

$ cd Replication
$ make -C cpp_tests

The cpp_tests directory contains a Makefile that can be edited by the user that compiles the C++ script. The Makefile inside directory cpp_tests accesses the user-defined environment variable $GUROBI_HOME. If it is set, then the code for the model with Covariates inside laspated.h is accessible. Otherwise, the script will not run the experiments that use covariates.

After the C++ code is successfully compiled, you can run the replication script with:

$ python replication_script.py

This command will generate the rectangular, hexagonal and by district discretizations using the LASPATED python module, run the experiments with the C++ code and process the results, generating the figures and tables presented in the Quality Control and Reuse Potential sections of this paper.

All results presented in the paper are saved in the directory replication_results, with subdirectories plots and tables. To access the raw data generated by C++, go to the directory cpp_tests/results, with subdirectories for each experiment made. Subdirectories ex1 and ex2 contain the results for the first and second model without covariates (Tests from sections 7.1.1 and 7.1.2). Subdirectory ex3 contains the results for the model with covariates (Table 2) and real_data contains the results for the experiments with real data from the Reuse Potential section.

7.1 Quality Control Tests with Artificial Data without Covariates

In these tests, points arrive in 2-dimensional space and over time according to a periodic non-homogeneous Poisson process. There is only one type of point, thus $|\mathcal{C}|=1$, and hence the type notation c is omitted. The region under consideration is $\mathcal{S}={[0,10]}^2$, as shown in Figure 11.

The region $\mathcal{S}$ is partitioned into two subsets $\mathcal{S}=\mathcal{B}\cup \mathcal{R}$, with $\mathcal{B}=[0,5]\times (5,10]\cup (5,10]\times [0,5]$ and $\mathcal{R}=[0,5]\times [0,5]\cup (5,10]\times (5,10]$.

Figure 11

7.1.1 Test 1

In this test, the rate function $\lambda :\mathcal{S}\times [0,T]\mapsto {\mathbb{R}}_{+}$ is different on $\mathcal{B}$ and on $\mathcal{R}$, and is periodic with period 2, as follows:

$\lambda (\mathrm{s},\mathrm{t})=\{\begin{array}{ccc}0.1\hfill & \text{if}& \mathrm{s}\in \mathcal{B} \text{and} \mathrm{t}\in (2k,2k+1] \text{for some} k\in \mathbb{N},\hfill \\ 0.5\hfill & \text{if}& \mathrm{s}\in \mathcal{B} \text{and} \mathrm{t}\in (2k–1,2k] \text{for some} k\in \mathbb{N},\hfill \\ 0.5\hfill & \text{if}& \mathrm{s}\in \mathcal{R} \text{and} \mathrm{t}\in (2k,2k+1] \text{for some} k\in \mathbb{N},\hfill \\ 0.1\hfill & \text{if}& \mathrm{s}\in \mathcal{R} \text{and} \mathrm{t}\in (2k–1,2k] \text{for some} k\in \mathbb{N}.\hfill \end{array}$

This rate function is represented in Figure 12(a) (for time intervals that start at odd times) and Figure 12(b) (for time intervals that start at even times).

Figure 12

The user knows the region $\mathcal{S}$ and that $|\mathcal{C}|=1$, but does not know about the subregions $\mathcal{B}$ and $\mathcal{R}$ that affect the intensity function $\lambda$, and does not know that $\lambda$ is periodic with period 2. For estimation purposes, the user discretizes $\mathcal{S}$ into 100 square zones of unit area each, as shown in Figure 11. Thus $\mathcal{I}={\mathcal{I}}_{\mathcal{B}}\cup {\mathcal{I}}_{\mathcal{R}}$ indexes two types of zones, but this is not known by the user:

• ${\mathcal{I}}_{\mathcal{B}}$ indexes 50 blue zones (with blue circles in their centers in Figure 11); there are 25 blue zones in the bottom right and 25 blue zones in the upper left of the region;

• ${\mathcal{I}}_{\mathcal{R}}$ indexes 50 red zones (with red squares in their centers in Figure 11); there are 25 red zones on the bottom left and 25 red zones in the upper right of the region.

Also, for estimation purposes, the user discretizes time into 28 time intervals of length 1 each. Thus, the estimates are denoted with ${\lambda }_{i,t}$ for $i\in \mathcal{I}$ and $t\in \mathcal{T}:=\{1,\mathrm{\dots },28\}$. Given arrival data M_i,t,n for $i\in \mathcal{I}$, $t\in \mathcal{T}$, and $n=1,\mathrm{\dots },N_{i,t}$, the regularized loss function in (6) is used to estimate ${\lambda }_{i,t}$. The penalty coefficients are $w_{\mathit{\text{i,t}}}=w>0$ when i,j are neighboring zones, and $w_{\mathit{\text{i,t}}}=0$ otherwise, where w will be varied as described later. Two zones are neighbors if their borders share an edge.

We consider two partitions ${\mathcal{G}}_2$ and ${\mathcal{G}}_4$ of time intervals $\mathcal{T}$ in (6). For partition ${\mathcal{G}}_2$, the time groups are $G_0=\{2k:k=1,\mathrm{\dots },14\}$ and $G_1=\{2k+1:k=0,1,\mathrm{\dots },13\}$. For partition ${\mathcal{G}}_4$, the time groups are $G_0=\{4k:k=1,\mathrm{\dots },7\}$, $G_1=\{4k+1:k=0,1,\mathrm{\dots },6\}$, $G_2=\{4k+2:k=0,1,\mathrm{\dots },6\}$, and $G_3=\{4k+3:k=0,1,\mathrm{\dots },6\}$. For both partitions, penalty coefficients $W_G=w$ for all groups G (note that all groups have the same cardinality so that it seems reasonable to choose the same weights for different groups). Also, $\epsilon =0.001$ in (17).

For any value of the penalty parameter w, let ${\widehat{\lambda }}_{i,t}^w$ denote the estimator of ${\lambda }_{i,t}$ produced by (8). Note that when the penalty parameter $w=0$, then ${\widehat{\lambda }}_{i,t}^0$ reduces to the empirical estimator of the intensities, which is the mean rate of arrivals in zone i and time interval t. For a range of values of the penalty parameter w, we computed the estimates ${\widehat{\lambda }}_{i,t}^w$, and Figure 13 shows the mean (over all zones i and time intervals t) relative error given by

$\frac{1}{|\mathcal{I}||\mathcal{T}|}\sum _{i\in \mathcal{I}}\sum _{t\in \mathcal{T}}\left|\frac{{\lambda }_{i,t}–{\widehat{\lambda }}_{i,t}^w}{{\lambda }_{i,t}}\right|$

as a function of penalty parameter w, for different values of N_i,t.Figure 13 shows the mean relative errors for the following four estimators: estimator with partition ${\mathcal{G}}_4$ of time intervals and with neighbor-based spatial regularization (legend “4 groups, neighbors” in the figure), estimator with partition ${\mathcal{G}}_4$ of time intervals and without spatial regularization (legend “4 groups, no neighbors” in the figure), estimator with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization (legend “2 groups, neighbors” in the figure), and estimator with partition ${\mathcal{G}}_2$ of time intervals and without spatial regularization (legend “2 groups, no neighbors” in the figure). Figure 13 presents results for four values of the sample size: $N_{i,t}=1$, $N_{i,t}=10$, $N_{i,t}=50$, and $N_{i,t}=500$. Figure 13 also shows, for the three larger sample sizes, the mean relative error obtained by choosing the penalty parameter w by cross validation as follows: The data are partitioned into 5 subsets. Then, for each of 5 replications we used one of the data subsets (a different subset for each replication, with 20% of data used for training) to compute ${\widehat{\lambda }}_{i,t}^w$ for a range of values of w, and we used the remaining data to compute the out-of-sample likelihood value for each value of w. Then we determined the penalty parameter $w^{\ast }$ that maximizes the average out-of-sample likelihood value over the 5 replications. Figure 13 shows the mean relative error for the resulting cross validation estimator with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization using the penalty parameter $w^{\ast }$. Note that the empirical estimator, with $w=0$, does not depend on the partition of time intervals and the neighbor structure, and therefore is given by all four estimators described above at $w=0$. Table 1 shows the minimum mean relative errors over different values of penalty parameter w, as well as the values of w that attain the minimum, for each of the four estimators, plus the mean relative errors of the cross validation estimator (using the penalty parameter $w^{\ast }$) and the empirical estimator ($w=0$), for each of the four sample sizes considered in Figure 13.

Figure 13

Table 1

The column headings denote the estimators as follows: Reg 1: estimator with partition ${\mathcal{G}}_4$ of time intervals and with neighbor-based spatial regularization; Reg 2: estimator with partition ${\mathcal{G}}_4$ of time intervals and without spatial regularization; Reg 3: estimator with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization; Reg 4: estimator with partition ${\mathcal{G}}_2$ of time intervals and without spatial regularization; CV: cross validation; Emp: empirical estimator. For each sample size and for each estimator Reg 1, Reg 2, Reg 3, Reg 4, each cell of the table gives two numbers separated by a slash (/): the first number is the minimum mean relative error over different values of penalty parameter w, and the second number is the value of w that attains the minimum. For the cross validation estimator, the first number is the mean relative error using the penalty parameter $w^{\ast }$, and the second number is the value of $w^{\ast }$. For the empirical estimator, the first number is the mean relative error using the penalty parameter $w=0$.

SAMPLE SIZE Ni,t	REG 1	REG 2	REG 3	REG 4	CV	EMP
1	0.47/1.2	0.67/2.0	0.34/0.6	0.51/1.4	-	1.539/0
10	0.31/0.06	0.30/0.2	0.26/0.06	0.17/0.2	0.42/0.2	0.543/0
50	0.18/0.008	0.12/0.04	0.15/0.004	0.08/0.036	0.37/0.04	0.256/0
500	0.08/0.0	0.03/0.004	0.08/0.0	0.02/0.0028	0.26/0.002	0.07687/0

In Figure 13, note that the vertical and horizontal scales of the sub-figures are different. As expected, the estimation error is decreasing in the sample size N_i,t. Note that the mean relative errors shown in Figure 13, and the values of the penalty parameter w with the smallest mean relative errors (shown with diamonds), were based on using the correct values ${\lambda }_{i,t}$ in the calculations, which are not known by the user. Even when regularization is not based on the correct subregions and partition of time intervals, it can result in better estimates than the empirical estimates (corresponding to $w=0$), as long as reasonable values are chosen for w. The results indicate that the improvement of the estimates using regularization is more pronounced when the number of observations is small: for $N_{i,t}=1$, the estimation error decreases from about 152% (for the empirical estimator) to about 35% for the best regularized estimator, and for $N_{i,t}=10$, the estimation error decreases from about 52% to about 19%. Therefore, when a small amount of data are available, the regularized estimator results in a smaller estimation error, even if the regularization is not based on the correct model structure, as long as relevant information and reasonable weights are used for the regularization.

Figure 14 shows the true intensities $\lambda (s,t)$ and the estimates ${\widehat{\lambda }}_{i,t}^{w^{\ast }}$ obtained with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization using the cross validation penalty parameters $w^{\ast }$, for $s\in \mathcal{R}$, zone $1\in {\mathcal{I}}_{\mathcal{R}}$, $N_{i,t}=1$ (top left plot), for $s\in \mathcal{R}$, zone $1\in {\mathcal{I}}_{\mathcal{R}}$, $N_{i,t}=10$ (bottom left plot), for $s\in \mathcal{B}$, zone $6\in {\mathcal{I}}_{\mathcal{B}}$, $N_{i,t}=1$ (top right plot), and for $s\in \mathcal{B}$, zone $6\in {\mathcal{I}}_{\mathcal{B}}$, $N_{i,t}=10$ (bottom right plot). It can be seen that if little data are available, then the regularized estimates are closer to the true intensities than the empirical estimates.

Figure 14

7.1.2 Test 2

In this test, we consider the same region $\mathcal{S}$ and the same subregions $\mathcal{B}$ and $\mathcal{R}$ as in Test 1. The intensity function $\lambda$ is also periodic, but not piecewise constant, as follows:

$\begin{array}{cc}& \lambda (x,y,t)\hfill \\ & =\{\begin{array}{ccc}\mathrm{x}+\mathrm{y}\hfill & \text{if}& (\mathrm{x},\mathrm{y})\in \mathcal{B} \text{and} \mathrm{t}\in (2k,2k+1] \text{for some} k\in \mathbb{N},\hfill \\ 5(\mathrm{x}+\mathrm{y})\hfill & \text{if}& (\mathrm{x},\mathrm{y})\in \mathcal{B} \text{and} \mathrm{t}\in (2k–1,2k] \text{for some} k\in \mathbb{N},\hfill \\ 5(\mathrm{x}+\mathrm{y})\hfill & \text{if}& (\mathrm{x},\mathrm{y})\in \mathcal{R} \text{and} \mathrm{t}\in (2k,2k+1] \text{for some} k\in \mathbb{N},\hfill \\ \mathrm{x}+\mathrm{y}\hfill & \text{if}& (\mathrm{x},\mathrm{y})\in \mathcal{R} \text{and} \mathrm{t}\in (2k–1,2k] \text{for some} k\in \mathbb{N}.\hfill \end{array}\hfill \end{array}$

This rate function is represented in Figure 15(a) (for time intervals that start at odd times) and Figure 15(b) (for time intervals that start at even times).

Figure 15

As in Test 1, the user knows the region $\mathcal{S}$ and that $|\mathcal{C}|=1$, but does not know about the subregions $\mathcal{B}$ and $\mathcal{R}$ that affect the intensity function $\lambda$, and does not know that $\lambda$ is periodic with period 2. For estimation purposes, the user discretizes $\mathcal{S}$ into 100 square zones of unit area each. The user computes the estimators ${\widehat{\lambda }}_{i,t}^w$ for each $i\in \mathcal{I}$ and each $t\in \mathcal{T}$, by solving optimization problem (16) using regularization with penalty parameter w and $\epsilon =0.001$. Each estimator ${\widehat{\lambda }}_{i,t}^w$ now approximates the rate ${\int }_{R_i}\lambda (x,y,t)\text{d}(x,y)$, where R_i denotes the square of unit area for the zone i. For example, for a zone i with $R_i=(a_i–1,a_i]\times (b_j–1,b_j]$ and $\lambda (x,y,t)=x+y$ for $(x,y)\in R_i$, it holds that

$\begin{array}{cc}\hfill {\int }_{R_i}\lambda (x,y,t)\text{d}(x,y)& ={\int }_{a_i–1}^{a_i}{\int }_{b_j–1}^{b_j}(x+y)\text{d}y\text{d}x=a_i+b_j–1\hfill \\ & =\lambda (x_i,y_i,t)=x_i+y_i\hfill \end{array}$

where $(x_i,y_i)=(a_i–0.5,b_j–0.5)$ is the centroid of the zone i.

The mean relative error of estimator ${\widehat{\lambda }}_{i,t}^w$ is now given by

$\begin{array}{cc}& \frac{1}{|\mathcal{I}||\mathcal{T}|}\sum _{i\in \mathcal{I}}\sum _{t\in \mathcal{T}}\frac{\left|{\displaystyle {\int }_{R_i}}\lambda (x,y,t)\text{d}(x,y)–{\widehat{\lambda }}_{i,t}^w\right|}{{\displaystyle {\int }_{R_i}}\lambda (x,y,t)\text{d}(x,y)}\hfill \\ & =\frac{1}{|\mathcal{I}||\mathcal{T}|}\sum _{i\in \mathcal{I}}\sum _{t\in \mathcal{T}}\frac{\left|\lambda (x_i,y_i,t)–{\widehat{\lambda }}_{i,t}^w\right|}{\lambda (x_i,y_i,t)}.\hfill \end{array}$

Parameter estimates ${\widehat{\lambda }}_{i,t}^w$ were computed for a range of values of the penalty parameter w, for the same four estimators as in Test 1: estimator with partition ${\mathcal{G}}_4$ of time intervals and with neighbor-based spatial regularization (legend “4 groups, neighbors” in Figure 16), estimator with partition ${\mathcal{G}}_4$ of time intervals and without spatial regularization (legend “4 groups, no neighbors” in Figure 16), estimator with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization (legend “2 groups, neighbors” in Figure 16), and estimator with partition ${\mathcal{G}}_2$ of time intervals and without spatial regularization (legend “2 groups, no neighbors” in Figure 16). As before, four values were used for the sample size: $N_{i,t}=1$, $N_{i,t}=10$, $N_{i,t}=50$, and $N_{i,t}=500$.Figure 16 shows the mean relative error for these four estimators, as well as the estimator with partition ${\mathcal{G}}_2$ of time intervals and with neighbor-based spatial regularization using the penalty parameter $w^{\ast }$ chosen with cross validation. As before, the empirical estimator corresponds to $w=0$. In this test intensities are very different between neighboring zones, and as a result space regularization is not helpful, but if little data are available then time regularization still provides better estimates than the empirical estimator.

Figure 16