# A python tool to determine the thickness of the hydrate layer around clinker
grains using SEM-BSE images.

Florian Kleiner  Franz Becker  Christiane Rößler  Horst-Michael Ludwig  F. A.
Finger-Institute for Building Materials Science, Bauhaus-University Weimar,
Germany  Faculty of Science - Chair of Mineralogy, University of Erlangen-
Nuernberg (FAU), Germany

###### Abstract

To accurately simulate the hydration process of cementitious materials,
understanding the growth rate of C-S-H layers around clinker grains is
crucial. Nonetheless, the thickness of the hydrate layer shows substantial
variation around individual grains, depending on their surrounding.
Consequently, it is not feasible to measure hydrate layers manually in a
reliable and reproducible manner. To address this challenge, a software has
been developed to statistically determine the C-S-H thickness, requiring
minimal manual interventions for thresholding and for setting limits like
particle size or circularity.

This study presents a tool, which automatically identifies suitable clinker
grains and and perform statistical measurements of their hydrate layer up to a
specimen age of 28 days. The findings reveal a significant increase in the
C-S-H layer, starting from 0.45 µm after 1 day and reaching 3.04 µm after 28
days. However, for older specimens, the measurement of the C-S-H layer was not
feasible due to limited pore space and clinker grains.

###### keywords:

hydration , calcium silicate hydrate , hydrate layer

††journal: Journal of Microscopy

##  1 Introduction

Understanding the formation of the micro-structure of cementitious materials
allows to model and thus improve their properties. Calcium-silicate-hydrate
(C-S-H) and portlandite (CH) are the main hydration products of ordinary
Portland cement. C-S-H is mainly responsible for the concrete strength. Hence,
its growth rate within the first 28 days is of major interest to model the
hydration process [1, 2].

The growth process can be observed in embedded and polished binder samples by
using a scanning electron microscope (SEM) in backscatter electron (BSE) mode
(see Figure 1). Analyzing the phase distribution by segmenting each
identifiable phase is a common practice [3]. Furthermore, it is fairly easy to
measure the particle size distribution of the alite (impure form of tricalcium
silicate, C3S) particles [4].

During the hydration of cementitious materials, C-S-H forms a layer around
clinker grains, and its thickness increases over time. This layer can be
differentiated into ($i$) the dense, inner product and ($ii$) the less dense,
needle like shaped outer product [5]. Especially the development of the outer
product is often limited by the surrounding particles or by other developing
phases like C-S-H or portlandite (CH). Bazzoni [6] demonstrated that is is
possible to determine the C-S-H needle length using images from a scanning
transmission electron microscope (STEM). Nevertheless, it is challenging to
measure the hydrate layer thickness (HLT) manually in a reliable and
reproducible manner. Hence, in this study, we established a workflow and
software tool to automatically determine the HLT.

##  2 Materials and Methods

In this study, a commercially available alite (MIII polymorph, Vustah, Czech
Republic) was utilized. The chemical composition (Table 1) was determined by
X-ray fluorescence spectroscopy (XRF), and the phase purity was confirmed at
99.4 ± 0.4 wt.-% (alite) using X-ray diffraction (XRD) analysis.

Table 1: Chemical composition of the utilized alite obtained by XRF.  Oxides | Alite in wt.-%
---|---
SiO2 | 26.34
TiO2 | 0.031
Al2O3 | 0.238
Fe2O3 | 0.094
Mn2O3 | 0.006
MgO | 1.700
CaO | 70.79
Na2O | 0.083

Alite pastes were prepared with a water-to-solid ratio of 0.5 and stored in
sealed cylindrical containers with a diameter of approx. 8 mm. The phase
composition and alite purity of similarly prepared samples stored for 1, 7,
and 28 days were assessed by quantitative XRD analysis using the G-factor
method [7]. The measurements were carried out using a Bruker D8 DaVinci
diffractometer (MA, USA).

After 1, 7, 14, 28, 84, and 365 days the hydration of the cylindrical prisms
was stopped by immersion in isopropanol and subsequent drying at 60∘C. The
specimens were cut, embedded in low-viscosity resin and mechanically polished
using an oil-based diamond paste with a grain size down to 0.25 µm. Finally,
the polished specimens were carbon coated (approx. 10 nm) to avoid charging
effects in the SEM. The specimens were transferred to a SEM (Helios G4 UX,
ThermoFisher Scientific, MA, USA) to capture high-resolution BSE images of
their surfaces using the built-in concentric backscatter detector (CBS). These
image sets were stitched to larger images. Afterwards, they were processed
using custom code, which is the main subject of this study. The code was
developed for Python 3.11 and was integrated into a Jupyter Notebook and it is
publicly available [8].

The image contrast was logarithmically enhanced using the library scikit-image
[9]. Subsequently, the images were denoised employing the Non Local Means
algorithm [10] implemented in the library OpenCV [11]. The raw and
preprocessed images have varying contrasts as shown in Figure 1A and B. Even
after preprocessing, the segmentation cannot be carried out using fixed
threshold values. Therefore, the images were segmented using two manually
selected threshold values into three phases: ($i$) pores, ($ii$) hydration
products (CH and C-S-H), and ($iii$) unhydrated clinker. The hydration
products were not separated since this is not always possible without major
segmentation errors. However, the software is also able to process pre-
segmented images (e.g. segmented by machine learning algorithms), if data is
provided in a supported format (3-channel TIFF, one colour channel per phase).
Finally, the holes within particles were removed to avoid any errors caused by
pores within clinker grains.

Figure 1: Detail of the stitched, unprocessed BSE images of embedded and
polished alite sections after 7 (A) and 28 days (B) of hydration. Identifiable
phases from dark to light grey: resin/pores, C-S-H, CH, alite. The images also
show a typical deviation in image contrast, even after using identical imaging
parameters (12 kV, 0.8 nA).

For the herein presented data, the segmented images were used to identify
clinker particles with a grain diameter between 0.3 and 9.0 µm (grain area
between 0.07 to 63.6 µm2). Furthermore, particles with a circularity lower
than 0.3 were excluded to avoid extreme sectioning effects. For very oval
particles, the hydrate fringe could otherwise be greatly overestimated and
undercuts could also lead to undesirable results.

In addition, the clinker particles shrink over time and small grains are
consumed. Since the selected grain diameter was kept constant regardless of
the specimen age, this effect is not taken into account in this analysis.

For each particle, a circular region with a radius of 7.0 µm was selected
(Figure 2A, B). This image was then transformed into polar coordinates (Figure
2C). Particles located too close to the image border (less than 7 µm) were
excluded from further processing. The HLT was determined using 360 steps,
representing a 1∘ angle per step. To achieve this, the distance from the alite
phase (left-hand side of Figure 2C, blue) to the next pore (red) was measured.
If there was no adjacent pore within the 7.0 µm radius, the measurement was
stopped (Figure 2D, areas with adjacent grey). To reduce errors, 3∘ before and
after those areas were ignored as well (Figure 2D, marked as darker vertical
lines).

Figure 2: Illustration of the processing that is done to measure the HLT. In
the center of A, B and E is an unhydrated clinker (1) grain within alite
hydrated for 7 days. A: raw image; B: segmented image (red: pores, green:
hydrates, blue: alite); C: transformation of B into polar coordinates, the
clinker grain (blue) is on the left; D: like C, but only with the identified
hydrate in green (grey areas were excluded from processing); E: version of D
transformed into Cartesian coordinates.

Subsequently, the histogram of the HLT distribution was analyzed. To safely
determine the maxima $l_{\text{max}}$, a function was fit to the HLT histogram
using equation (1), where $l$ is the hydrate layer thickness and $f$ is the
frequency.

$f=a\cdot e^{-0.5\left(\frac{log(l-d)-b}{c}\right)}$ (1)

Additionally, by using the area of the alite $A_{\text{C}_{3}\text{S}}$ and
the hydrate phases $A_{\text{hydrate}}$, obtained from the image segmentation,
it was possible to estimate the amount of hydrate phase $\sigma_{\text{hyd}}$
(CH + C-S-H) using Equation (2). The densities used to calculate these values
were 3.14 g/cm3 ($\rho_{\text{C}_{3}\text{S}}$), 2.2 g/cm3 (CH) and $\approx$
2.0 g/cm3 (C-S-H) [5, 12]. For the hydrate mixture, a volume ratio C-S-H:CH of
2.5:1 was assumed, giving an approximate hydrate density $\rho_{\text{hyd}}$
of 2.07 g/cm3.

$\sigma_{\text{hydrate}}=\frac{A_{\text{hydrate}}\cdot\rho_{\text{hydrate}}}{A_{\text{hydrate}}\cdot\rho_{\text{hydrate}}+A_{\text{C}_{3}\text{S}}\cdot\rho_{\text{C}_{3}\text{S}}}$
(2)

##  3 Results and Discussion

The herein presented method was applied to alite specimens of different ages
in order to test its applicability. However, an existing dataset (1 day) and
datasets explicitly created for this research (7, 14, 28, 84, and 365 days)
were mixed. The datasets differ mainly in the imaged area and slightly in SEM
imaging parameters. For example, the accelerating voltage used varied from 10
to 15 kV. An overview of the images used is shown in Table 2. The raw images
are publicly available [13]

Table 2: Basic information and selected results of the alite datasets processed. Image size, normed particle count (diameter < 9.0 µm, circularity > 0.3), phase areas based on the segmentation, fit parameters ($a$-$d$), and the maximum HLT $l_{\text{max}}$ of the second value accumulation.  age in days | 1 | 7 | 14 | 28 | 84 | 365
---|---|---|---|---|---|---
area in mm2 | 0.56 | 2.41 | 2.41 | 2.41 | 2.41 | 2.41
particles |  |  |  |  |  |
per mm2 | 25296 | 6540 | 17719 | 8050 | 8385 | 7512
pores in |  |  |  |  |  |
area-% | 38.5 | 28.5 | 22.2 | 13.7 | 21.7 | 14.8
hydrates in |  |  |  |  |  |
area-% | 38.0 | 52.8 | 44.5 | 78.7 | 72.7 | 83.2
clinker in |  |  |  |  |  |
area-% | 23.9 | 19.1 | 33.6 | 7.7 | 5.7 | 2.0
$a$ in µm | 0.999 | 0.929 | 0.898 | 0.898 | - | -
$b$ |  $-$1.506 | 0.097 | 0.585 | 0.858 | - | -
$c$ | 1.003 | 0.471 | 0.487 | 0.475 | - | -
$d$ | 0.219 | 0.297 | 0.680 |  $-$1.000 | - | -
$l_{\text{max}}$ in µm | 0.45 | 1.43 | 1.69 | 3.04 | - | -

On close inspection of the images, it becomes evident that the low-viscosity
resin did not completely penetrate the specimens. This is apparent in Figure
2A (2), where the large pore on the top center appears a little darker than
the surface of the epoxy resin in other parts of the image. In the most cases,
this is not an issue, as those holes are interpreted as pores anyway. However,
larger unfilled pores should be manually painted black, since this could lead
to erroneously identified clinker particles or hydrates. Additionally, the
needle-like structure of the C-S-H and the penetration depth of the electron
beam (from 1.60 µm in alite up to 3.75 µm in epoxy resin at 12 kV [14]) cause
a brightness gradient between pores and hydrate phase, particularly in
denoised images.

A significant challenge lies in accurately determining the threshold value
that distinguishes the pore space from the hydrate phase (see (3) in Figure 2B
and C). The polishing process leaves the specimen surface with slight
unevenness due to differences in material hardness, resulting in undesirable
topography and brightness variations at phase borders. This segmentation
artifact, combined with cracks within the clinker particles might cause the
first peak of the HLT measurements (see orange distributions in Figure 3).
Figure 2B (4) also indicates that in some instances, edges of hydrate were
incorrectly identified as clinker phase, which may also affect the later
analysis. These issues are further enhanced, because the datasets were
stitched together from several smaller images, whose brightness may show
slight deviations.

If the HLT is compared to the grain diameters (Figure 3), it is noticeable
that the data density decreases with increasing age for the same examined area
size. This is indicated by an increasing level of noise, particularly in the
case of larger grain diameters. This phenomenon can be attributed to the
decrease in the amount of unhydrated particles present in older specimens (see
Table 2), leading to shifts in the particle size distribution over time. The
grey histogram above the 2D-distribution, which shows the logarithmic
frequency of measurements per grain diameter bin, also indicates this
alteration.

These 2D-distributions indicate accumulations of HLT measurements, depending
on the alite grain diameters. To eliminate the afore mentioned effect of the
particle size distribution, the the maximum amount of hydrate layer
measurements was normalized to 1 for each particle diameter step. These
normalized distributions were then used to calculate the cumulative HLT
distribution depicted in orange on the right of the 2D-distributions in Figure
3.

Figure 3: Normalized distributions of the HLT of hydrated alite (aged 1, 7,
14, 28, 84, and 365 days) depending on the particle diameter. The histogram on
top illustrates the logarithmic frequency for each diameter bin in the diagram
below. The histogram on the right illustrates the frequency of HLT
measurements as sum of the main diagram. The values are thus independent of
the frequency of single grain diameters.

All distributions show an accumulation at very low measurements (< 0.5 to 1.0
µm), which could be an artifact caused by the preparation and segmentation.
For example, a very small hydrate fringe appears to be around some particles
after the automatic segmentation (compare Figure 2B c, (3)). However, this
layer would not exist after a manual segmentation.

Histograms of specimens with an age up to 28 days show two clusters of HLTs.
However, only the second accumulation seems to be a valid result. Therefore,
the first accumulation was ignored. For older specimens this can be observed
only for larger particles (e.g. a diameter above 2.0 µm for 14 and 28 days
hydrated alite), while smaller particles only seem to have the first
accumulation. Nevertheless, a trend of an increasing HLT between 1 and 28 days
can be observed.

To obtain the value of $l_{\text{max}}$, Equation (1) was fitted to the orange
HLT histogram in Figure 3. The resulting curves are plotted in red. However,
due to the absence of a second accumulation in specimens older than 28 days,
fitting was not feasible in those cases. The resulting fit parameters $a$,
$b$, $c$, $d$, and $l_{\text{max}}$ are listed in Table 2.

The selected limit value for circularity of 0.3 was a compromise between a
high data density for larger particles and signal clarity. Increasing this
value can improve the signal clarity of the second accumulation, if there are
enough suitable particles available. Furthermore, this reduces the chance of
overestimating the HLT due to angular distortions for non-circular particles.
However, this also reduces the amount of measures, especially for particles
with a diameter larger than 5.0 µm. This is an issue for older specimens,
where the particle count is significantly reduced due to the hydration
progress. For example, in order to be able to make any statement at all about
specimens hydrated for 84 and 365 days, the circularity value must be lowered
to at least 0.4 to obtain enough data for particles larger than 5.0 µm in
diameter. Nevertheless, they do not seem to have a second cluster of HLTs in
the range measured.

There are multiple possible reasons for this result, which can be followed on
the basis of Figure 4. Firstly, the pore space between hydrates is
significantly smaller, reducing the possibility to obtain good measurements.
Secondly, the few remaining large particles appear to disintegrate into
several smaller particles or break up due to sample preparation, leaving some
apparent pores. Finally, a large part of the remaining space is filled by CH.
Therefore, in this case the results should be considered invalid.

Figure 4: Image detail of 365 days hydrated alite. This image shows a lack of
pore space filled by CH, disintegrated larger particles and significantly
smaller particles than in Figure 2 (note the scale bar differences).

It is important not to confuse the increase of the HLT with the hydration
progress. Figure 5 shows the hydration progress, measured using XRD (dashed
lines). The combined ratio of C-S-H and CH is represented by the green dash-
dotted line, the HLT is plotted in a solid red line. While both, the XRD and
HLT measurements indicate an increase, the HLT does not follow the same
pattern as the XRD measurements.

Figure 5: Comparisson between the hydration progression measured using XRD
(dashed lines) with $l_{\text{max}}$ (solid line). The pink markers (X)
indicate the amount of hydrate phase $\sigma_{\text{hyd}}$ estimated from the
image segmentation.

This can be argued with multiple reasons:

  * 1.

First and foremost, completely dissolved grains are not taken into account to
measure the hydrate layer thickness. This concerns a significant area of the
image and becomes more prevalent with specimen age.

  * 2.

The grains could be cut at an unfavourable position (far off-center), which
leads to an apparently larger hydrate layer.

  * 3.

CH and C-S-H phases are not separated. However, this is not anticipated to
pose a significant problem, as the majority of measurements will be
automatically excluded due to the absence or minimal presence of pores within
the 7.0 µm radius, particularly when the particle is embedded in CH.

In comparison, Masoero et al. [15] calculated a reaction zone of about 0.4 to
0.5 µm. They used different particle size distributions of alite in a
suspension (water/solid ratio of 50). In their proposed model, the reaction
zone should be filled after 21 h. This seems to be in the same range of the
measurement of about 0.45 µm in this study.

In contrast, for sieved C3S (mean particle size 6 µm) Costoya Fernández [16]
calculated a peak thickness of C-S-H of 0.3 µm after 12 h of hydration.
Powders with larger particle size distributions showed a significantly smaller
C-S-H thickness (C-S-H thickness of 0.1 µm for a mean particle size of 18 µm).

Garrault et al. [5] proposed a maximum dense C-S-H-layer thickness of 0.4 µm,
with a tendency to even lower values for larger grain diameters. Bazzoni [6]
measured the C-S-H needles on top of the dense C-S-H and found, they were in
the range of 0.3 to 0.4 µm in length after 7 days. Combining both results, the
overall HLT should therefore not exceed 0.8 µm, which is not supported by the
data found in this study.

In some studies [15, 16, 5], the C-S-H layer thickness was calculated from ion
concentrations and heat flow measurements. They are therefore based on
volumetric measurements of the whole specimens and not individual alite
grains. In contrast, our study only includes C-S-H layers directly around a
particle and ignores a bulk of the formed hydrate.

The results of the estimated hydrate phase based on image segmentation are
shown in Figure 5 (pink markers (X)). They do not give a clear pattern, but
usually the amount of hydrate was overestimated. This could be due to the fact
that $\sigma_{\text{hydrate}}$ was calculated from areas $A_{\text{hydrate}}$
and $A_{\text{C}_{3}\text{S}}$ instead of their respective volumes and since
only a rough phase ratio and density estimate was used. However, this can not
explain the deviation in the data for the 14 day specimen. Therefore, it is
likely that not enough area was measured for a reliable phase measurement.

Finally, the computational requirements and processing time of the script used
are non-negligible factors. The images with an area of 2.4 µm2 were scaled
from a pixel size of 0.0422 µm to 0.0843 µm to reduce processing time and
memory requirements. The processing of the resized image requires about 8 GB
of RAM and can take up to 12 hours on an Intel Core i7 8700, depending on the
image size and the amount of particles to be processed. However, there is
still some optimization potential by improving the multi-threaded processing
or by compiling the code to machine language.

##  4 Conclusions and outlook

The herein presented method can be utilized to determine the development of
the HLT of alite for young to medium aged specimens (up to 28 days) in an
objective manner. Aside from the image segmentation itself, it is a fully
automated process, requiring little user intervention.

The C-S-H growth with increasing age is clearly reflected in the results. The
maxima in the HLT distribution increases from 0.45 µm after 1 day up and
reaching 3.04 µm after 28 days. However, at 84 days or later the HLT
distribution does not allow to evaluate the HLT distribution using the
selected parameters. Consequently, the proposed approach finds applicability
when the microstructure is not overly dense. In practical terms, this
signifies the method’s viability requires the presence of a substantial amount
of pores and clinker particles for its successful application.

Additionally, it was tried to calculate the hydration progress from the
generated data. However, BSE images up to 2.4 mm2 in size still seem to be too
small for a reliable evaluation.

While the script was only applied to a simplified cementitious system (clean
alite). First tests not presented in this paper did indicate that this method
may also be applicable for clean belite specimens. Furthermore, it should be
expandable to more complex multiphase materials. However, this imposes higher
demands on the segmentation quality. Applying more advanced segmentation
methods (e.g. machine learning based) could therefore improve the reliability
of the results.

In the future, the influence of parameters like the circularity limit should
be studied using synthetic datasets. Furthermore, multiple images of specimens
of varying sizes but the same age should be studied, to determine the
specimen-to-specimen variability. Finally, this method could be improved to
determine the HLT of 3-dimensional, voxel-based datasets.

## Author contributions

Florian Kleiner: Conceptualization, Methodology, Investigation, Software,
Visualization, Data curation, Writing – original draft. Franz Becker:
Supporting measurements (XRF, XRD), Writing – review & editing, Validation
Christiane Rößler: Supervision, Writing – review & editing, Validation. Horst-
Michael Ludwig: Financing, Supervision, Writing – review & editing,
Validation.

## Acknowledgment

The research was supported by the Deutsche Forschungsgemeinschaft (DFG), grant
number 344069666. Open access funding enabled and organized by project DEAL.
Thanks to M. Böhme, S. Unbehau and F. Ellermann for their helpful input.

## Conflict of interest

The authors declare no potential conflict of interests.

## Supporting information

Additional supporting information may be found in the online version of the
article at the publisher’s website.

## References

  * [1] L. Nguyen-Tuan, C. Rößler, H.-M. Ludwig, Simulation of C3S Hydration at the Micro-scale using Level Set Method, in: 41st Cement & Concrete Science Conference, 2022.
  * [2] L. Nguyen-Tuan, M. A. Etzold, C. Rößler, H.-M. Ludwig, Growth and porosity of C-S-H phases using the sheet growth model, Cement and Concrete Research 129 (2020) 105960.  doi:10.1016/j.cemconres.2019.105960.
  * [3] D. C. Chu, J. Kleib, M. Amar, M. Benzerzour, N.-E. Abriak, Determination of the degree of hydration of Portland cement using three different approaches: Scanning electron microscopy (SEM-BSE) and Thermogravimetric analysis (TGA), Case Studies in Construction Materials 15 (2021) e00754.  doi:10.1016/j.cscm.2021.e00754.
  * [4] K. Lyu, W. She, C. Miao, H. Chang, Y. Gu, Quantitative characterization of pore morphology in hardened cement paste via SEM-BSE image analysis, Construction and Building Materials 202 (2019) 589–602.  doi:10.1016/j.conbuildmat.2019.01.055.
  * [5] S. Garrault, T. Behr, A. Nonat, Formation of the C-S-H Layer during early hydration of tricalcium silicate grains with different sizes, The journal of physical chemistry. B 110 (1) (2006) 270–275.  doi:10.1021/jp0547212.
  * [6] A. Bazzoni, Study of early hydration mechanisms of cement by means of electron microscopy (2014).  doi:10.5075/EPFL-THESIS-6296.
  * [7] D. Jansen, F. Goetz-Neunhoeffer, C. Stabler, J. Neubauer, A remastered external standard method applied to the quantification of early OPC hydration, Cement and Concrete Research 41 (6) (2011) 602–608.  doi:10.1016/j.cemconres.2011.03.004.
  * [8] F. Kleiner, kleinerELM/measure_hydrate_layer (2023).  doi:10.5281/ZENODO.8211756.
  * [9] S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, T. Yu, scikit-image: image processing in Python, PeerJ 2 (2014) e453.  doi:10.7717/peerj.453.
  * [10] A. Buades, B. Coll, J.-M. Morel, Non-Local Means Denoising, Image Processing On Line 1 (2011) 208–212.  doi:10.5201/ipol.2011.bcm_nlm.
  * [11] G. Bradski, The OpenCV Library, Dr. Dobb’s Journal of Software Tools (2000).
  * [12] H. M. Jennings, A model for the microstructure of calcium silicate hydrate in cement paste, Cement and Concrete Research 30 (1) (2000) 101–116.  doi:10.1016/S0008-8846(99)00209-4.
  * [13] F. Kleiner, Large SEM-BSE images of hydrated alite of ages from 1 day up to 1 year (2023).  doi:10.5281/ZENODO.8210797.
  * [14] F. Kleiner, C. Matthes, C. Rößler, Argon broad ion beam sectioning and high resolution scanning electron microscopy imaging of hydrated alite, Cement and Concrete Research 150 (2021) 106583.  doi:10.1016/j.cemconres.2021.106583.
  * [15] E. Masoero, J. J. Thomas, H. M. Jennings, A Reaction Zone Hypothesis for the Effects of Particle Size and Water-to-Cement Ratio on the Early Hydration Kinetics of C 3 S, Journal of the American Ceramic Society 97 (3) (2014) 967–975.  doi:10.1111/jace.12713.
  * [16] M. M. Costoya Fernández, Effect of particle size on the hydration kinetics and microstructural development of tricalcium silicate, Ph.D. thesis (2008).  doi:10.5075/EPFL-THESIS-4102.


# Addressing Bias in Active Learning with Depth Uncertainty Networks… or Not

Chelsea Murray James U. Allingham Javier Antorán José Miguel Hernández-Lobato
Department of Engineering
University of Cambridge
{clm88, jua23, ja666<EMAIL_ADDRESS>
###### Abstract

Farquhar et al. (2021) show that correcting for active learning bias with
underparameterised models leads to improved downstream performance. For
overparameterised models such as NNs, however, correction leads either to
decreased or unchanged performance. They suggest that this is due to an
“overfitting bias” which offsets the active learning bias. We show that depth
uncertainty networks operate in a low overfitting regime, much like
underparameterised models. They should therefore see an increase in
performance with bias correction. Surprisingly, they do not. We propose that
this negative result, as well as the results Farquhar et al. (2021), can be
explained via the lens of the bias-variance decomposition of generalisation
error.

##  1 Introduction

Active learning improves data efficiency by identifying which examples
maximise the expected gain in model performance (Cohn et al., 1995; Settles,
2010). An acknowledged but relatively poorly understood problem with this
approach, however, is that actively selecting informative training points
introduces bias in inference, as the training data no longer follow the
population distribution (MacKay, 1992). Given that many statistical results
rely on data points being identically and independently distributed (i.i.d.)
samples from the population distribution, applying standard estimators to
actively sampled datasets results in optimising for the wrong objective
(Farquhar et al., 2021).

Recently, Farquhar et al. (2021) introduced an estimator,
$\tilde{R}_{\text{LURE}}$, that corrects for the bias induced by active
sampling by re-weighting the contribution of individual data points to the
estimated population risk. Farquhar et al. (2021) find that, while
underparameterised models such as linear regression benefit from the
elimination of active learning bias, the performance of overparameterised
models such as neural networks is either unaffected or negatively impacted by
the unbiased risk estimator. This result is explained as active learning bias
acting as a regulariser and helping to prevent flexible models from
overfitting on small datasets.

Depth uncertainty networks (DUNs) are a class of Bayesian neural network (BNN)
in which a prior is placed over the depth of the network, allowing us to
favour simpler (i.e., shallower) functions. This enables DUNs to adapt their
complexity to the varying size of the training dataset as additional labels
are acquired during active learning. Because of this, we posit that DUNs will
suffer less from overfitting in the early stages of training than other BNNs.
This would imply that using the unbiased risk estimator of Farquhar et al.
(2021) to eliminate active sampling bias should improve the performance of
DUNs. Our main contributions are:

  * •

Empirically validating the hypothesis that DUNs present less overfitting bias
than other BNNs, specifically Monte Carlo dropout (MCDO), during active
learning;

  * •

Investigating the application of the $\tilde{R}_{\text{LURE}}$ estimator to
eliminate active learning bias from DUNs. Unexpectedly, despite DUNs having
larger active learning bias than overfitting bias, we find that
$\tilde{R}_{\text{LURE}}$ does not result in improved downstream performance;

  * •

Providing hypotheses for why we do not observe an impact on DUNs’ downstream
performance when applying the unbiased risk estimator.

##  2 Background

###  2.1 Bias in active learning

We use notation similar to Farquhar et al. (2021). Given a loss function
$\mathcal{L}(\mathbf{y},f_{\theta}(\mathbf{x}))$, we aim to find $\theta$ that
minimises the population risk over $p_{\text{data}}(\mathbf{y},\mathbf{x})$:
${\mathbb{E}_{\mathbf{x},\mathbf{y}\sim p_{\text{data
}}}\left[\mathcal{L}\left(\mathbf{y},f_{\theta}(\mathbf{x})\right)\right]}$.
In practice, we only have access to $N$ samples from $p_{\text{data}}$. These
yield the empirical risk,
${r=\frac{1}{N}\sum_{n=1}^{N}\mathcal{L}\left(\mathbf{y}_{n},f_{\theta}\left(\mathbf{x}_{n}\right)\right)}$.
$r$ is an unbiased estimator of the population risk when the data are drawn
i.i.d. from $p_{\text{data}}$. If these samples have been used to train
$\theta$, our estimators become biased. We refer to estimators biased by
computing with the training data with capital letters (i.e., $R$). In the
active learning setting, our model is optimised using a subset of $M$ actively
sampled data points:

$\tilde{R}=\frac{1}{M}\sum_{m=1}^{M}\mathcal{L}\left(\mathbf{y}_{i_{m}},f_{\theta}\left(\mathbf{x}_{i_{m}}\right)\right),\quad
i_{m}\sim\alpha\left(i_{1:m-1},\mathcal{D}_{\text{pool }}\right).$ (1)

The proposal distribution $\alpha\left(i_{m};i_{1:m-1},\mathcal{D}_{\text{pool
}}\right)$ over the pool of unlabelled data $\mathcal{D}_{\text{pool}}$
represents the probability of each index being sampled next, given that we
have already acquired
$\mathcal{D}_{\text{train}}=\\{\mathbf{x}_{i}\\}_{1:m-1}$. As a result of
active sampling, $\tilde{R}$ is not an unbiased estimator of $R$. Farquhar et
al. (2021) propose $\tilde{R}_{\text{LURE}}$ (“levelled unbiased risk
estimator”), which removes active learning bias:

$\displaystyle\tilde{R}_{\text{LURE}}\\!\equiv\\!\frac{1}{M}\sum_{m=1}^{M}v_{m}\mathcal{L}_{i_{m}};\quad
v_{m}\\!\equiv 1+\frac{N-M}{N-m}\left(\frac{1}{(N-m+1)\
\alpha\left(i_{m};i_{1:m-1},\mathcal{D}_{\text{pool}}\right)}-1\right),$ (2)

where
$\mathcal{L}_{i_{m}}\equiv\mathcal{L}\left(\mathbf{y}_{i_{m}},f_{\theta}\left(\mathbf{x}_{i_{m}}\right)\right)$.
Intuitively, the estimator works by re-weighting each example’s contribution
to the total loss by its inverse acquisition probability, such that if unusual
examples have particularly high acquisition probabilities, they contribute
proportionally less to the total risk. A detailed derivation of this estimator
is given by Farquhar et al. (2021).

Surprisingly, Farquhar et al. (2021) find that using $\tilde{R}_{\text{LURE}}$
during training can negatively impact the performance of flexible models. They
explain this by noting that the overfitting bias ($r{-}R$) typical in such
models has the opposite sign to the bias induced by active learning. That is,
active learning encourages the selection of “difficult” to predict data
points, increasing the overall risk, whereas overfitting results in the total
risk being underestimated. Where overfitting is substantial, it is (partially)
offset by active learning bias, such that correcting for active learning bias
does not improve model performance.

###  2.2 Depth uncertainty networks for Active Learning

$\mathbf{x}^{(n)}$$\mathbf{y}^{(n)}$$\theta\ $$d$$\beta$$N$ Figure 1: DUN
graphical model.

DUNs, depicted in figure 1, are a BNN variant in which the depth of the
network $d\in[1,2,\\!...,\\!D]$ is treated as a random variable. Model weights
$\theta$ are kept deterministic, simplifying optimisation. We place a
categorical prior over depth
$p(d)\,{=}\,\mathrm{Cat}(d|\\{\beta_{i}\\}_{i=0}^{D})$, resulting in a
categorical posterior which can be marginalised in closed form. This procedure
can be seen as Bayesian model averaging over an ensemble of subnetworks of
increasing depth. Thus, for DUNs, the predictive distribution and Bayesian
active learning by disagreement (BALD) acquisition objective (Houlsby et al.,
2011) are tractable and cheap to compute. See appendix A or (Antorán et al.,
2020) for a detailed description.

In DUNs, predictions from different depths of the network correspond to
different kinds of functions—shallow networks induce simple functions, while
deeper networks induce more complex functions. We choose a prior that favours
shallow models, limiting the flexibility of DUNs when few datapoints are
observed. As the actively acquired set increases, the model likelihood will
start to dominate the prior, leading the posterior over depth to favour deeper
models and thus more flexible functions. We show this in Figure C.6. We
hypothesize that the capability to automatically adapt model flexibility as
more data are observed will reduce overfitting in the small data regimes
typical of active learning. Therefore DUNs, unlike standard BNNs, may in fact
benefit from the use active learning bias corrective weights. It is worth
noting that automatically adapting model complexity to the observed data is
also a promised feature of traditional BNNs. However, unlike in DUNs where
inferences are exact, weight space models require crude approximations which
fall short in this regard (Ghosh et al., 2019).

##  3 Results

First, using the analysis of Farquhar et al. (2021), we show that the weighted
estimator removes active learning bias for DUNs. We then show that DUNs
overfit less than MCDO. Finally, we examine how training with
$\tilde{R}_{\text{LURE}}$ impacts downstream performance for DUNs. Experiments
are performed on nine UCI regression datasets (Hernández-Lobato and Adams,
2015); results for the Boston, Concrete and Energy datasets are presented in
the following sections, with results for the rest provided in appendix C. The
experimental setup is described in appendix B. All experiments are repeated 40
times, with the mean and standard deviations over the repetitions reported in
our figures.

###  3.1 Quantifying active learning bias

We estimate the extent of active learning bias when using $\tilde{r}$ and
verify the unbiasedness of $\tilde{r}_{\text{LURE}}$. A DUN is trained on
$1,000$ randomly selected data points from $\mathcal{D}_{\text{pool}}$ using
the standard negative log-likelihood (NLL) loss. We then estimate the model’s
risk using both $\tilde{r}$ and $\tilde{r}_{\text{LURE}}$, drawing $M$ samples
from an unobserved test set $\mathcal{D}_{\text{test}}$. The active learning
bias is estimated as

$B_{\text{ALB}}(\bullet)=r-E_{x_{1:M}\sim\alpha(x_{1:M};\mathcal{D}_{\text{pool}})}[\bullet],$
(3)

where $\bullet$ is either $\tilde{r}$ or $\tilde{r}_{\text{LURE}}$ and $r$ is
computed on the whole $\mathcal{D}_{\text{test}}$. Figure 2 shows that
$\tilde{r}_{\text{LURE}}$ is unbiased, whereas the bias in $\tilde{r}$
approaches zero only as $M\to N$.

Figure 2: Active learning bias in risk estimation with DUNs, evaluated using
$\tilde{r}$ (blue) and $\tilde{r}_{\text{LURE}}$ (orange). The unbiased risk
estimator $\tilde{r}_{\text{LURE}}$ eliminates active learning bias, while the
standard estimator $\tilde{r}$ tends to overestimate risk due to $\alpha$
prioritising more difficult to fit points.

###  3.2 Overfitting bias

To test the hypothesis that DUNs’ prior over depth reduces overfitting in low
data regimes, we follow the procedure of Farquhar et al. (2021) to isolate
overfitting bias. We train a model with $M$ actively sampled points using
$\tilde{R}_{\text{LURE}}$ as the training objective. After acquiring each new
point, we evaluate

$B_{\text{OFB}}=r-\tilde{R}_{\text{LURE}}.$ (4)

The difference between the risk evaluated on the test set ($r$) and on the
train set once the active learning bias is removed ($\tilde{R}_{\text{LURE}}$)
is taken as a measure of overfitting ($B_{\text{OFB}}$).

Figure 3 shows that the magnitude of $B_{\text{OFB}}$ for DUNs is small
relative to BNNs trained with MCDO (Gal and Ghahramani, 2016), one of the
methods used in Farquhar et al. (2021). This suggests that DUNs are, indeed,
robust to overfitting on small datasets. Overfitting bias in DUNs is also
larger in magnitude than active learning bias for the datasets shown
(comparing figure 3 to figure 2). This is not the case for MCDO (comparing
figure 3 to figure C.2).

Figure 3: Evolution of overfitting bias for DUNs (orange) and MCDO (green)
trained with $\tilde{R}_{\text{LURE}}$ as the amount of acquired data $M$
grows. Overfitting bias is smaller in DUNs than for MCDO.

###  3.3 Training with unbiased risk estimators

According to Farquhar et al. (2021), removing the active learning bias in
cases where ${|B_{\text{ALB}}|\gg|B_{\text{OFB}}|}$ should positively affect
performance. Indeed, they observe this for linear models. Thus, having
confirmed that the LURE weights eliminate active learning bias in DUNs, and
that DUNs are in the $|{B_{\text{ALB}}|\gg|B_{\text{OFB}}}|$ regime, at least
for certain datasets, we evaluate whether using the unbiased estimator during
training improves downstream performance. We compare the NLL obtained by
evaluating our models on $\mathcal{D}_{\text{test}}$ when using $\tilde{R}$
and $\tilde{R}_{\text{LURE}}$ as the training objectives. Surprisingly, figure
4 shows that using $\tilde{R}_{\text{LURE}}$ does not affect NLL performance
at all.

Figure 4: Test NLL as more training points are acquired for DUNs trained with
$\tilde{R}$ (blue) and $\tilde{R}_{\text{LURE}}$ (orange). Active learning
bias correction has minimal impact on downstream performance.

##  4 Discussion

In figure 4, we measure our models’ generalisation error, given by
$\mathtt{bias}^{2}+\mathtt{variance}$ (Bishop, 2007). Considering that neural
networks are very flexible, we expect model bias to be negligible, and active
learning bias to dominate. When we eliminate the latter with importance
weights, as in equation 2, we pay the price of additional variance (Owen,
2013). We therefore hypothesise that our results can be explained by gains due
to bias reduction being lost due to increased variance. This increase in
variance will depend on the choice of sampling distribution $\alpha$.
“Overfitting bias”, as presented in equation 4 and (Farquhar et al., 2021) is
computed using the train set and thus is not directly comparable with active
learning bias in equation 3. Furthermore, the effect of “overfitting bias” is
captured by model variance. The error induced by model variance is always
additive with any statistical biases that may exist, precluding any
cancellation effects. This hypothesis also explains why Farquhar et al. (2021)
find $\tilde{R}_{\text{LURE}}$ to be successful for linear models: since
inflexible models are well-specified by small amounts of data, applying
$\tilde{R}_{\text{LURE}}$ barely increases variance. We also highlight that in
the experimental design used by Farquhar et al. (2021) to obtain their linear
model results, the temperature of the proposal distribution ($T=10,000$) is
much larger than that used here ($T=10$). Higher temperatures approach a
deterministic proposal, interfering with the assumption that $\alpha(\cdot)$
has full support.

## Acknowledgments

JA acknowledges support from Microsoft Research, through its PhD Scholarship
Programme, and from the EPSRC. JUA acknowledges funding from the EPSRC and the
Michael E. Fisher Studentship in Machine Learning. CM acknowledges assistance
from Sebastian Farquhar in clarifying the implementation of the active
learning bias experiments. This work has been performed using resources
provided by the Cambridge Tier-2 system operated by the University of
Cambridge Research Computing Service (http://www.hpc.cam.ac.uk) funded by
EPSRC Tier-2 capital grant EP/P020259/1.

## References

  * Antorán et al. [2020] Javier Antorán, James Allingham, and José Miguel Hernández-Lobato.  Depth uncertainty in neural networks.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Bishop [2007] Christopher M. Bishop.  _Pattern Recognition and Machine Learning, 5th Edition_.  Information science and statistics. Springer, 2007.
  * Cohn et al. [1995] David A Cohn, Zoubin Ghahramani, and Michael I Jordan.  Active learning with statistical models.  Technical report, MASSACHUSETTS INST OF TECH CAMBRIDGE ARTIFICIAL INTELLIGENCE LAB, 1995.
  * D’Amour et al. [2020] A D’Amour, K Heller, D Moldovan, B Adlam, B Alipanahi, A Beutel, C Chen, J Deaton, J Eisenstein, MD Hoffman, et al.  Underspecification presents challenges for credibility in modern machine learning. arxiv 2020.  _arXiv preprint arXiv:2011.03395_ , 2020.
  * Farquhar et al. [2021] Sebastian Farquhar, Yarin Gal, and Tom Rainforth.  On statistical bias in active learning: How and when to fix it.  _International Conference on Learning Representations_ , 2021.
  * Gal and Ghahramani [2016] Yarin Gal and Zoubin Ghahramani.  Dropout as a Bayesian approximation: Representing model uncertainty in deep learning.  In _International Conference on Machine Learning_ , pages 1050–1059. PMLR, 2016.
  * Ghosh et al. [2019] Soumya Ghosh, Jiayu Yao, and Finale Doshi-Velez.  Model selection in bayesian neural networks via horseshoe priors.  _Journal of Machine Learning Research_ , 20(182):1–46, 2019.  URL http://jmlr.org/papers/v20/19-236.html.
  * Hernández-Lobato and Adams [2015] José Miguel Hernández-Lobato and Ryan Adams.  Probabilistic backpropagation for scalable learning of Bayesian neural networks.  In _International Conference on Machine Learning_ , pages 1861–1869. PMLR, 2015.
  * Houlsby et al. [2011] Neil Houlsby, Ferenc Huszar, Zoubin Ghahramani, and Máté Lengyel.  Bayesian active learning for classification and preference learning.  _CoRR_ , abs/1112.5745, 2011.  URL http://arxiv.org/abs/1112.5745.
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _arXiv preprint arXiv:1502.03167_ , 2015.
  * MacKay [1992] David JC MacKay.  Information-based objective functions for active data selection.  _Neural computation_ , 4(4):590–604, 1992.
  * Owen [2013] Art B. Owen.  _Monte Carlo theory, methods and examples_.  2013\.
  * Settles [2010] Burr Settles.  Active learning literature survey.  _Machine Learning_ , 15(2):201–221, 2010.
  * Zeiler and Fergus [2014] Matthew D. Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In David J. Fleet, Tomás Pajdla, Bernt Schiele, and Tinne Tuytelaars, editors, _Computer Vision - ECCV 2014 - 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part I_ , volume 8689 of _Lecture Notes in Computer Science_ , pages 818–833. Springer, 2014.

##  Appendix A Depth uncertainty networks

BNNs translate uncertainty in weight-space into predictive uncertainty by
marginalising out the posterior over weights. The intractability of this
posterior necessitates the use of approximate inference methods. DUNs, in
contrast, leverage uncertainty about the appropriate depth of the network in
order to obtain predictive uncertainty estimates. There are two primary
advantages of this approach: 1) the posterior over depth is tractable,
mitigating the need for limiting approximations; and 2) due to the sequential
structure of feed-forward neural networks, both inference and prediction can
be performed with a single forward pass, making DUNs suitable for deployment
in low-resource environments [Antorán et al., 2020].

DUNs are composed of subnetworks of increasing depth, with each subnetwork
contributing one prediction to the final model. The computational model for
DUNs is shown in figure A.5. Inputs are passed through an input block,
$f_{0}(\cdot)$, and then sequentially through each of $D$ intermediate blocks
$\\{f_{i}(\cdot)\\}_{i=1}^{D}$, with the activations of the previous layer
acting as the inputs of the current layer. The activations of each of the
$D+1$ layers are passed through the output block $f_{D+1}(\cdot)$ to generate
per-depth predictions. These are combined via marginalisation over the depth
posterior, equivalent to forming an ensemble with networks of increasing
depth. Variation between the predictions from each depth can be used to obtain
predictive uncertainty estimates, in much the same way that variance in the
predictions of different sampled weight configurations yield predictive
uncertainties in BNNs.

$f_{0}$ $f_{1}$ $f_{2}$ $f_{3}$ $f_{4}$ $f_{D}$ $\mathbf{x}$ $f_{D+1}$
$\mathbf{\hat{y}}_{i}$ Figure A.5: DUN computational model [Antorán et al.,
2020]. Each layer’s activations are passed through the output block, producing
per-depth predictions.

###  A.1 Uncertainty over depth

Figure A.6 compares the graphical model representations of a BNN and a DUN. In
BNNs, the weights $\theta$ are treated as random variables drawn from a
distribution $p_{\gamma}(\theta)$ with hyperparameters $\gamma$. In DUNs, the
depth of the network $d$ is assumed to be stochastic, while the weights are
learned as hyperparameters. A categorical prior distribution is assumed for
$d$, with hyperparameters $\beta$:
$p_{\beta}(d)=\operatorname{Cat}(d\mid\\{\beta_{i}\\}_{i=0}^{D})$. The model’s
marginal log likelihood (MLL) is given by

$\log
p(\mathcal{D}_{\text{train}};\theta)=\log\sum_{i=0}^{D}\left(p_{\beta}(d=i)\prod_{n=1}^{N}p\left(\mathbf{y}_{n}\mid\mathbf{x}_{n},d=i;\theta\right)\right),$
(5)

where the likelihood for each depth is parameterised by the corresponding
subnetwork’s output:
$p(\mathbf{y}\mid\mathbf{x},d=i;\theta)=p\left(\mathbf{y}\mid
f_{D+1}\left(\mathbf{a}_{i};\theta\right)\right)$, where
$\mathbf{a}_{i}=f_{i}(\mathbf{a}_{i-1})$.

$\mathbf{x}^{(n)}$$\mathbf{y}^{(n)}$$\theta$$\gamma$$d$$N$

$\mathbf{x}^{(n)}$$\mathbf{y}^{(n)}$$\theta\ $$d$$\beta$$N$

Figure A.6: Left: graphical model of a BNN. Right: graphical model of a DUN.

The posterior,

$p(d\mid\mathcal{D}_{\text{train}};\theta)=\frac{p(\mathcal{D}_{\text{train}}\mid
d;\theta)p_{\beta}(d)}{p(\mathcal{D}_{\text{train}};\theta)},$ (6)

is also a categorical distribution, whose probabilities reflect how well each
depth subnetwork explains the data. DUNs leverage two properties of modern
neural networks: first, that successive layers have been shown to extract
features at increasing levels of abstraction [Zeiler and Fergus, 2014]; and
second, that networks are typically underspecified, meaning that several
different models may explain the data well [D’Amour et al., 2020]. The first
property implies that initial layers in a network specialise on low-level
feature extraction, yielding poor predictions from the shallower subnetworks
in a DUN. This is handled by the depth posterior assigning low probabilities
to earlier layers, in preference for later layers that specialise on
prediction. The second property suggests that subnetworks at several depths
are able to explain the data simultaneously and thus have non-zero posterior
probabilities, yielding the diversity in predictions required to obtain useful
estimates of model uncertainty.

###  A.2 Inference in DUNs

Antorán et al. [2020] find that directly optimising the MLL equation 5 with
respect to $\theta$ results in convergence to a local optimum in which all but
one layer is attributed a posterior probability of zero. That is, direct
optimisation returns a deterministic network of arbitrary depth. The authors
instead propose a stochastic gradient variational inference (VI) approach with
the aim of separating optimisation of the weights $\theta$ from the posterior.
A surrogate categorical distribution $q_{\phi}(d)$ is introduced as the
variational posterior. The following evidence lower bound (ELBO) can be
derived:

$\mathcal{L}(\phi,\theta)=\sum_{n=1}^{N}\mathbb{E}_{q_{\phi}(d)}\left[\log
p\left(\mathbf{y}_{n}\mid\mathbf{x}_{n},d;\theta\right)\right]-\mathrm{KL}\left(q_{\phi}(d)\|p_{\beta}(d)\right),$
(7)

allowing $\theta$ and the variational parameters $\phi$ to be optimised
simultaneously. Antorán et al. [2020] show that under the VI approach, the
variational parameters converge more slowly than the weights, enabling the
weights to reach a setting at which a positive posterior probability is learnt
for several depths.

It is important to note that equation 7 can be computed exactly, and that
optimising $\mathcal{L}(\phi,\theta)$ is equivalent to exact inference of the
true posterior $p(d\mid\mathcal{D}_{\text{train}};\theta)$ in the limit: since
both $q$ and $p$ are categorical, equation 7 is convex and tight at the
optimum (i.e., $q_{\phi}(d)=p(d\mid\mathcal{D}_{\text{train}};\theta)$). Since
the expectation term can be computed exactly using the activations at each
layer (recall $p(\mathbf{y}\mid\mathbf{x},d=i;\theta)=p\left(\mathbf{y}\mid
f_{D+1}\left(\mathbf{a}_{i};\theta\right)\right)$), the ELBO can be evaluated
exactly without Monte Carlo sampling.

##  Appendix B Experimental setup

###  B.1 Data

We implement batch-based active learning, with batches of $20$ data points
acquired in each query (with the exception of Yacht, which has a query size of
10 due to its smaller size). An initial training set representing a small
proportion of the full data is selected uniformly at random in the first
query. For all datasets, 80% of the data are used for training, 10% for
validation and 10% for testing. Details about dataset sizes, input
dimensionality, and active learning specifications for each dataset are
provided in table B.1.

Table B.1: Summary of datasets and active learning specifications. 80% of the
data are used for training, 10% for validation and 10% for testing.

Name | Size | Input Dim. | Init. train size | No. queries | Query size
---|---|---|---|---|---
Boston Housing | 506 | 13 | 20 | 17 | 20
Concrete Strength | 1,030 | 8 | 50 | 30 | 20
Energy Efficiency | 768 | 8 | 50 | 30 | 20
Kin8nm | 8,192 | 8 | 50 | 30 | 20
Naval Propulsion | 11,934 | 16 | 50 | 30 | 20
Power Plant | 9,568 | 4 | 50 | 30 | 20
Protein Structure | 45,730 | 9 | 50 | 30 | 20
Wine Quality Red | 1,599 | 11 | 50 | 30 | 20
Yacht Hydrodynamics | 308 | 6 | 20 | 20 | 10

###  B.2 Model specifications

We implement a fully-connected network with residual connections, with 100
hidden nodes per layer. The networks contain 10 hidden layers for DUNs, or
three hidden layers for MCDO. We use ReLU activations, and for DUNs batch
normalisation is applied after every layer [Ioffe and Szegedy, 2015].
Optimisation is performed over $1,000$ iterations with full-batch gradient
descent, with momentum of $0.9$ and a learning rate of $10^{-4}$. A weight
decay value of $10^{-5}$ is also used. We do not implement early stopping, but
the best model based on evaluation of the evidence lower bound on the
validation set is used for reporting final results. For MCDO models a fixed
dropout probability of $0.1$ is used and prediction is based on 10 MC samples.
DUNs use a uniform categorical prior over depth. The hyperparameter settings
largely follow those used in Antorán et al. [2020] for their experiments with
toy regression datasets.

All experiments are repeated 40 times with different weight initialisations
and train-test splits (with the exception of the Protein dataset, for which
experiments are repeated only 30 times due to the cost of evaluation on the
larger test set). Unless otherwise specified, we report the mean and standard
deviation of the relevant metric over the repeated experiment runs. Following
Farquhar et al. [2021], the proposal distribution
$\alpha\left(i_{m};i_{1:m-1},\mathcal{D}_{\text{pool }}\right)$ used is a
stochastic relaxation of the Bayesian active learning by disagreement (BALD)
acquisition function proposed by [Houlsby et al., 2011]. BALD scores are
multiplied by a constant temperature $T$ before being passed through the
`softmax` function; candidate points are then sampled with probability
proportional to the `softmax` probabilities. We use $T=10$ as this is found to
yield the best performance for most of the datasets considered (results for
experiments to determine the optimal value of $T$ are given in section C.4).

##  Appendix C Additional experimental results

###  C.1 Quantifying active learning bias

Figure C.1 presents the results of the active learning bias experiments, as in
figure 2, for the remaining datasets. Figure C.2 shows equivalent results for
BNNs trained with MCDO.

Figure C.1: Bias introduced by actively sampling points with DUNs, evaluated
using $\tilde{R}$ (blue) and $\tilde{R}_{\text{LURE}}$ (orange).

Figure C.2: Bias introduced by actively sampling points with MCDO, evaluated
using $\tilde{R}$ (blue) and $\tilde{R}_{\text{LURE}}$ (orange).

###  C.2 Overfitting bias

Figure C.3 compares the magnitude of overfitting bias for DUNs and BNNs
trained with MCDO, as in figure 3. Figure C.3 shows additionally the
overfitting bias when training with $\tilde{R}$.

Figure C.3: Overfitting bias for DUNs (solid lines) and MCDO (dashed lines)
trained with $\tilde{R}$ (blue) and $\tilde{R}_{\text{LURE}}$ (orange). NLL
for Protein and Yacht is displayed in log scale.

###  C.3 Training with unbiased risk estimators

Figure C.4 presents the results of training DUNs with $\tilde{R}$ and
$\tilde{R}_{\text{LURE}}$, as in figure 4, for the remaining datasets.

Figure C.4: Test NLL for DUNs trained with $\tilde{R}$ (blue) and
$\tilde{R}_{\text{LURE}}$ (orange). A larger value of the orange line
indicates that training with $\tilde{R}_{\text{LURE}}$ harms performance,
despite removing active learning bias.

###  C.4 Temperature of the proposal distribution

Figure C.5 shows the test NLL for DUNs using different temperatures $T$ for
the proposal distribution.

Figure C.5: Test NLL of DUNs using a stochastic relaxation of BALD. Different
temperatures of the proposal distribution are compared. Means of 40 runs of
the experiments are shown; standard deviations are not plotted for clarity.

###  C.5 Depth posteriors for DUNs at different stages of active learning

Figure C.6 shows the posterior distributions over depth for DUNs for the
smallest and largest dataset sizes used during active learning.

Figure C.6: Posterior probabilities over depth for DUNs trained on UCI
regression datasets, for the smallest (blue bars) and largest (orange bars)
labelled datasets used in active learning.


# Changes in Commuter Behavior from COVID-19 Lockdowns in the Atlanta
Metropolitan Area

Tejas Santanam H. Milton Stewart School of Industrial and Systems Engineering,
Georgia Institute of Technology Anthony Trasatti H. Milton Stewart School of
Industrial and Systems Engineering, Georgia Institute of Technology Hanyu
Zhang H. Milton Stewart School of Industrial and Systems Engineering, Georgia
Institute of Technology Connor Riley H. Milton Stewart School of Industrial
and Systems Engineering, Georgia Institute of Technology Pascal Van Hentenryck
H. Milton Stewart School of Industrial and Systems Engineering, Georgia
Institute of Technology Ramayya Krishnan Heinz College, Carnegie Mellon
University

###### Abstract

This paper analyzes the impact of COVID-19 related lockdowns in the Atlanta,
Georgia metropolitan area by examining commuter patterns in three periods:
prior to, during, and after the pandemic lockdown. A cellular phone location
dataset is utilized in a novel pipeline to infer the home and work locations
of thousands of users from the Density-based Spatial Clustering of
Applications with Noise (DBSCAN) algorithm. The coordinates derived from the
clustering are put through a reverse geocoding process from which word
embeddings are extracted in order to categorize the industry of each work
place based on the workplace name and Point of Interest (POI) mapping.
Frequencies of commute from home locations to work locations are analyzed in
and across all three time periods. Public health and economic factors are
discussed to explain potential reasons for the observed changes in commuter
patterns.

##  1 Introduction

This paper is motivated by a desire to gain an understanding of the impact of
COVID-19 on work-related travel behaviors in the Atlanta metropolitan area.
The paper analyzes and answers the following questions:

  1. 1.

Where do people in the Atlanta metropolitan area live and work?

  2. 2.

How did people’s commuter travel behaviors change over the course of the
COVID-19 pandemic and its associated lockdowns?

  3. 3.

How have different industries changed their work styles as a result of the
COVID-19 pandemic and its associated lockdowns?

From a high-level perspective, this paper provides a novel methodology for
answering the above questions that consists of the following steps:

  1. 1.

Divide cell phone location data to users active before, during, and after
COVID-19 lockdowns in the Atlanta area.

  2. 2.

Apply density-based clustering to identify home locations, work locations, and
days present at each.

  3. 3.

Use reverse geocoding combined with a bag-of-words model to map locations to
workplace addresses and industries.

  4. 4.

Analyze user commuter and work behaviors in relation to work industry before,
during, and after COVID-19 lockdowns in the Atlanta area.

The paper is organized as follows. Section 2 first presents a short survey of
related literature. Next, Section 3 describes the cell phone location dataset
and its characteristics. Section 4 follows with details on the clustering
procedure, while Section 5 delves into the reverse geocoding and address
mapping step. Section 6 presents the final results of the analysis and Section
7 ends with a discussion of the results.

##  2 Literature Review

There have been a few recent papers that use cell phone location, GPS data,
and smart card data to understand urban mobility patterns. The use of cell
phone location data to generate home and work trip pairs, as well as assign
some trip purposes are explored by [6]. [1] compare Estonian mobile data with
population register data to analyze the effectiveness of a model that
determined home and work locations for users. More recently, [20] use smart
card data for home and work location detection, while [24] use mobile data to
understand the spatial structure of urban commuting. [23] also use mobile data
to examine how people move around under the context of key activity locations.
Patterns in home-work commuting are also examined by [13]. [28] focus on
learning taxi driver behavior from GPS data.

The use of density-based clustering methods on mobile phone data has also been
explored in prior literature. [11] improve on the existing DBSCAN algorithm
methods to identify trip ends from Chinese mobile phone location data. [27]
use DBSCAN clustering to cluster public transit passengers into regional
travel areas. [18] implement DBSCAN for noise reduction as they infer route
choice for commuter trips. DBSCAN is used as the final step in deriving the
characteristics of commuting patterns on key commuting routes by [25].

A few recent papers have looked at the impact of the COVID-19 pandemic on
human mobility patterns. [7] estimate how individual mobility outside the home
is affected by local disease prevalence and restriction orders to stay-at-
home. [3] document the evolution of work-from-home behaviors over the course
of the pandemic and the expectation for that kind of working behavior in the
future. [4] uses cell phone location data to understand the impact of the
pandemic on nursing home entries. [10] look at the correlation between
demographic and socioeconomic variables and home-dwelling time records as a
result of stay-at-home orders derived from large-scale mobile phone location
tracking data. [26] build an interactive tool and process anonymized mobile
device location data to identify trips and produce a set of variables
including a social distancing index, the percentage of people staying at home,
the number of visits to work and non-work locations, out-of-town trips, and
trip distance. [12] look at how the number of COVID-19 cases can be used to
predict the following week’s travel behaviors in Houston, Texas. [17] analyze
the impacts of COVID-19 on micromobility use and perception.

The choice of the Atlanta metropolitan area as a case study in this paper is
also based in prior literature analyzing location data, commuter behavior, and
commuter location data. [21] use GPS data to analyze destination choice
behavior for vehicle trips in Atlanta. [19] and [22] discuss the consistency
of commuter behavior using public rail in Atlanta. The impact on this
consistency in Atlanta as a result of the pandemic is analyzed by [2] who also
propose a resilient multimodal transit system to serve the Atlanta area in a
pandemic case in the context of COVID-19. [9] consider commuter home locations
in the integration of latent demand into the design of an on-demand multimodal
transit system for Atlanta. The use of the OpenStreetMap Application
Programming Interface (API) for Atlanta area points of interest has also been
implemented by [14] who build a tool to explore the transit accessibility in
the Atlanta metropolitan area to those points of interest.

Based on the literature above, there exists a gap in using cell phone (or any
similar type of) location data to derive a breakdown of remote and in-person
commute behaviors by work industry. Ultimately, the main novel contributions
of this paper relate to the analysis of lockdown effects on commuter
behaviors. To our knowledge, there are no comparable methods in the existing
literature for deriving granular analysis on commuting patterns by industry
from mobile location data. The pipeline that connects density-based clustering
to reverse geocoding to work type is also a unique and novel contribution, and
is scalable in its potential applications outside of the Atlanta area and in
temporal periods beyond the focus on COVID-19 lockdowns.

##  3 Data

The paper uses a dataset from X-Mode which contains the cell phone latitude
and longitude location data for a set of anonymized users who use mobile
applications with the X-Mode platform integrated. X-Mode is a data collection
platform for location data that app developers can choose to integrate into
their app [15]. The original dataset contains over 12 terabytes of location
data for users all over the United States, so the first step is to preprocess
the data to focus on the target time and locations. The data is filtered and
contained to a bounding box around the metropolitan area surrounding the city
of Atlanta. The city of Atlanta is located in Fulton County, Georgia which is
the county in Georgia that saw the largest number of COVID-19 cases and
related deaths in 2020, with 56,800 reported COVID-19 cases over the year. The
resulting data points are further filtered to include user locations during
the months of January 2020 (baseline pre-pandemic), April 2020 (pandemic
lockdown), and July 2020 (pandemic post-lockdown). The Gross Domestic Product
(GDP) of the state of Georgia grew by 5.1% in January 2020 compared to January
2019, decreased by 16.6% in April 2020 compared to April 2019 (the largest
decrease in the last decade), and fell by 5.5% in July 2020 compared to July
2019. The monthly GDP trend is shown in Figure 1. These months are analyzed
and focused on in this paper.

Figure 1: Monthly change in Georgia GDP (early 2020)

In January 2020, which was used as a baseline, there were around 430,000
unique users who had at least one location point tracked and logged to the
dataset. In order to achieve a generalization of the users’ travel patterns,
the number of active days for each user in January is calculated. A threshold
of greater than 24 days with at least one data point present in the dataset is
set to maintain a strong clustering result. With this threshold, a little over
21,000 users are identified as having enough information in the month of
January. The location data for a single month for one anonymized user is
displayed on the map as point clusters and a heat map in Figure 2.

Figure 2: Location data for one user over a month

##  4 Clustering Methods

To begin the analysis of commuter behavior from the location data, work and
home locations are needed to be determined. To do so, a clustering approach is
used. Work location and home location for each user are identified using
DBSCAN clustering. DBSCAN is a density-based clustering algorithm that takes a
set of points in some space and groups together points that are closely packed
together [8]. It is used here due to its quality performance in finding core
samples of high density at home and work locations and expanding clusters from
them. The following process of location identification below is used:

  1. 1.

Starting with the filtered data containing 21,000 users from above, further
filter the data to the group of users with more than 300 data points of
activity during the pre-pandemic period in January. If a user has fewer than
300 activities, the home and work locations derived might not be accurate.
Those users with fewer than 300 data points of activity are not assigned any
home and work locations. 4,500 of the 21,000 user subset meet this criteria.

  2. 2.

Taking the filtered group of users, divide the user’s GPS data into workdays
and weekend. For each sub-dataset, cluster the GPS data using DBSCAN. The
parameters for DBSCAN are set to be eps = 0.001, and minsample = 5% of the
sample size and were manually tuned.

  3. 3.

Find the cluster centers for the above users for workdays and weekends.

  4. 4.

For each cluster center, find which days they appear in within a month.

  5. 5.

The cluster center that appears most frequently in terms of days in January
(that also has a frequency greater than 17 days) is classified as the “home
location”. If there is no “home location” found, then go to the next user.

  6. 6.

Remove the “home location” from the list of cluster centers and find the most
frequent location that appear in workdays in the list of cluster centers. If
that location appears more than 10 days per month, then identify that location
as “work location”.

Figure 3: Home and Work locations in the 1.5K user subset Figure 4: Home to
Work distances in the 1.5K user subset

Using this method, the anonymous user in Figure 2 will have their home
location identified within the top cluster of 1,658 points and their work
location identified within the leftmost cluster of 371 points. After repeating
the above steps on the 21,000 user subset, 1,428 users are also above the 24
day activity threshold in April and July as well. This roughly 1.5K user
subset is the main focus of the analysis going forward in this paper. These
clusters were also manually inspected to obtain a high-level of confidence in
the accuracy of the home and work labels based on the general locations of
where one might expect to see residential and commercial areas. The home
locations (in black) and work locations (in blue) of some of those 1,428 users
are shown as markers and marker clusters in Figure 3. Home and work locations
for any particular user are connected in the figure by a red line. The
distribution of the distances between the home and work locations for these
users are shown in kilometers in Figure 4.

##  5 Reverse Geocoding

To figure out the work types of these users and analyze the change in their
work patterns during the pandemic, a reverse geocoding procedure is applied on
the latitude and longitude of their work locations. The Nominatim
OpenStreetMap API is used to retrieve the actual name of the work locations;
for example, “Drury Inn” or “Wells Fargo” [16] [5]. For the set of work
locations where the Nominatim OpenStreetMap API is unable to provide a
workplace name, a database of Atlanta points of interest is used to map the
work locations to workplace names. For the remainder of the few unmapped work
locations, a manual mapping is done through Google Maps.

After the workplace names for all of the users are gathered, a bag-of-words
model is run on the work place names. The purpose of the bag-of-words model is
to count the frequency of the words that appear in the names of workplaces in
the roughly 1.5K workplace set. Some high-frequency words are identified such
as “school”, “hotel”, and “hospital”. These high-frequency words serve as
guidance for designing the word embeddings to extract categories from location
names. The generated categories are presented in Table 1. The word embeddings
extract categories from the location names. For example, any workplace with
the word “University”, “School”, or “Academy“ in the title is mapped by the
word embeddings program to the Education category (except in the case of
workplace names like “Emory University Hospital” which is mapped to the
Medical category). Similarly, workplace names containing “Taco”, “Pizza”, or
“Diner” are mapped to the Restaurant category.

Figure 5: Category of Nominatim OpenStreetMap API address matches

##  6 Results

Table 1 displays the number of people associated with each workplace category
out of the 1,428 users analyzed. The notes column details the types of
workplaces that were mapped to each category at the reverse geocoding step.
The largest categories contained people in Mix, Retail, White Collar, and
Education jobs. A comparison between the 1,428 user subset and the roughly
21,000 users who were only active in January in the dataset (not including the
users in the 1,428 user subset) is created, and it is found that the
distribution of job type by percentage is roughly similar. This can verify
that the subsets generated in the analysis are representative of the overall
Atlanta metropolitan area community. The median workdays over each month and
the associated percentiles are shown in Table 2. In that table, the 25%
percentiles refers to fewer than 25% of the users in that category working
fewer than the specified number of days. The median number of days is graphed
in Figure 6. While all occupation types see a decrease in the number of days
users went in to work, the Blue Collar and Medical fields see the greatest
resilience of in-person work. Some fields like the White Collar and Education
sectors see their median day count drop to zero as a result of the possibility
of virtual transition these sectors have. Some sectors like the hotel industry
or entertainment venues also see a median day count drop to zero as a result
of forced closings and layoffs from lockdown policies in Atlanta. Figure 7
shows a comparison between the percentage of days worked in person during the
pandemic months compared to the expected percentage of days worked in person
from the North American Industry Classification System (NAICS) dataset, which
is a government dataset for businesses in Canada, Mexico, and the United
States on how necessary in-person work is within different occupations and
industries.

Categories | Count | Note
---|---|---
Airport | 14 | Airport Terminals, Parking Lots
Blue Collar | 59 |  Manufacturing, Vehicle Repair, Cemetery Work
Education | 156 | K-12 schools and Higher education, Preschools
Entertainment | 9 |  Stadium, theatre, theater, movie, bowling, opera, museum, orchestra
Government | 30 | City Halls, Mints, Consulates
Hotel | 51 | Hotels, Motels
Medical | 130 |  Hospitals, Urgent Care, Dentists/Chiropractic Offices
Mix | 385 | Small Businesses, Home Businesses, etc.
Recreation | 16 | Trails, Parks, Golf Courses, Tennis Courts
Religion | 28 | Churches, Temples, Mosques
Residential | 70 | Apartments, Townhomes, Condos
Restaurant | 90 | Food and Drink eat-in Establishments
Retail | 203 |  Brick-and-mortar Retail
White Collar | 187 | Large Office Buildings/Companies

Table 1: Amount of users by workplace category in the 1.5K user subset

| Jan | Apr | July
---|---|---|---
Categories | 50% | 25% | 75% | 50% | 25% | 75% | 50% | 25% | 75%
Other | 20 | 17 | 22 | 4 | 0 | 17 | 7 | 0 | 19
Airport | 17 | 16 | 22 | 1 | 0 | 4 | 4 | 0 | 13
Blue Collar | 19 | 16 | 21 | 14 | 1 | 20 | 16 | 3 | 20
Education | 18 | 16 | 19 | 0 | 0 | 1 | 0 | 0 | 3
Entertainment | 18 | 16 | 19 | 0 | 0 | 0 | 1 | 0 | 4
Government | 18 | 15 | 19 | 1 | 0 | 3 | 1 | 0 | 5
Hotel | 20 | 16 | 22 | 0 | 0 | 7 | 0 | 0 | 15
Medical | 19 | 16 | 21 | 10 | 1 | 19 | 12 | 2 | 18
Recreation | 18 | 14 | 22 | 0 | 0 | 3 | 2 | 0 | 18
Religion | 20 | 15 | 23 | 1 | 0 | 6 | 7 | 0 | 13
Residential | 20 | 17 | 23 | 10 | 0 | 22 | 4 | 0 | 22
Restaurant | 21 | 18 | 24 | 5 | 0 | 18 | 15 | 1 | 20
Retail | 21 | 18 | 24 | 6 | 0 | 20 | 14 | 1 | 20
White Collar | 19 | 16 | 21 | 0 | 0 | 2 | 0 | 0 | 5

Table 2: Median and quartiles of workdays by workplace category Figure 6: Bar
graph of median workdays by workplace category Figure 7: Comparison of actual
workday percentage to NAICS in-person work percentages

##  7 Discussion

The most impactful conclusion stems from the three different types of median
workday changes over time seen in the results. One type of trajectory is a
sudden large drop in in-person work as the COVID-19 pandemic struck that
continues post-lockdown as well. This is seen in categories like White Collar,
Education, Entertainment, Hotel, and Government. The second type of trajectory
that is noticeable is a large drop as COVID-19 struck that then bounces back
(but not all the way to normal levels) post-lockdown. This trajectory is seen
in categories like Restaurants and Retail locations. The third type of
trajectory is a small drop to a level that stays consistent through lockdown
and post-lockdown, as seen in Blue Collar and Medical work. For the most part,
the first category contains jobs that either can easily be transitioned
virtually or jobs that are impossible to transition virtually and thus had to
pause. The second category contains those kinds of jobs that had to
temporarily stop during lockdown but were able to come back in non-traditional
ways such as outdoor seating-only for restaurants or online order and pick-up
operations in retail locations. The third category contains those jobs that
had to continue in-person regardless of the pandemic conditions, such as
medical work or running a grocery store. Further conclusions can be drawn from
the NAICS percentages. While some in-person work percentages align with the
government expectation for work-from-home pre-COVID-19, as seen in the case of
Blue Collar work, many occupations see a large drop from the expectation
during lockdown, such as a 73% drop in Restaurant work and Airport work or a
97% drop in Hotel work. This large drop from the expectation can be attributed
to the effect of legislation and policy on in-person work. Since legislation
forced the closing of many service-based companies such as bowling alleys,
theaters and restaurants that did not have outdoor seating, these jobs that
would normally be considered by the government to be in-person were no longer
possible. Potential areas of future work include expanding this analysis to
other cities and comparing the results with Atlanta. Another interesting line
of analysis would be to analyze these workplaces from the customer perspective
and see where people went during COVID-19. Credit card data could be used to
compare where consumers spent their money in-person versus online. The
addition of smart card data and other forms of public transit data sources can
also be used to further the granularity of the analysis over time in this
paper. After some time, analysis of the long-term effects and implications of
the pandemic on the greater commuter culture would be useful to explore.

## References

  * [1] Rein Ahas, Siiri Silm, Erki Saluveer and Olle Järv  “Modelling home and work locations of populations using passive mobile positioning data”  In _Location based services and TeleCartography II_ Springer, 2009, pp. 301–315
  * [2] Ramon Auad et al.  “Resiliency of on-demand multimodal transit systems during a pandemic”  In _Transportation Research Part C: Emerging Technologies_ 133 Elsevier, 2021, pp. 103418
  * [3] Alexander Bick, Adam Blandin and Karel Mertens  “Work from home before and after the Covid-19 outbreak”  In _Available at SSRN 3786142_ , 2021
  * [4] M Keith Chen, Judith A Chevalier and Elisa F Long  “Nursing home staff networks and COVID-19”  In _Proceedings of the National Academy of Sciences_ 118.1 National Acad Sciences, 2021, pp. e2015455118
  * [5] Konstantin Clemens  “Geocoding with openstreetmap data”  In _GEOProcessing 2015_ , 2015, pp. 10
  * [6] Serdar Çolak et al.  “Analyzing cell phone location data for urban travel: current methods, limitations, and opportunities”  In _Transportation Research Record_ 2526.1 SAGE Publications Sage CA: Los Angeles, CA, 2015, pp. 126–135
  * [7] Samuel Engle, John Stromme and Anson Zhou  “Staying at home: mobility effects of covid-19”  In _Available at SSRN 3565703_ , 2020
  * [8] Martin Ester, Hans-Peter Kriegel, Jörg Sander and Xiaowei Xu  “A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise”  In _Proc. of 2nd International Conference on Knowledge Discovery and_ , 1996, pp. 226–231
  * [9] Hongzhao Guan, Beste Basciftci and Pascal Van Hentenryck  “Heuristic Algorithms for Integrating Latent Demand into the Design of Large-Scale On-Demand Multimodal Transit Systems”  In _arXiv preprint arXiv:2212.03460_ , 2022
  * [10] Xiao Huang et al.  “Staying at home is a privilege: Evidence from fine-grained mobile phone location data in the United States during the COVID-19 pandemic”  In _Annals of the American Association of Geographers_ 112.1 Taylor & Francis, 2022, pp. 286–305
  * [11] Haihang Jiang et al.  “Improved F-DBSCAN for Trip End Identification Using Mobile Phone Data in Combination with Base Station Density”  In _Journal of Advanced Transportation_ 2022 Hindawi, 2022
  * [12] Junfeng Jiao, Mira Bhat and Amin Azimian  “Measuring travel behavior in Houston, Texas with mobility data during the 2020 COVID-19 outbreak”  In _Transportation letters_ 13.5-6 Taylor & Francis, 2021, pp. 461–472
  * [13] Kevin S Kung, Kael Greco, Stanislav Sobolevsky and Carlo Ratti  “Exploring universal patterns in human home-work commuting from mobile phone data”  In _PloS one_ 9.6 Public Library of Science San Francisco, USA, 2014, pp. e96180
  * [14] Alexander Li et al.  “Ease and Equity of Point of Interest Accessibility via Public Transit in the US”  In _arXiv preprint arXiv:2212.06954_ , 2022
  * [15] “Mode - empowering innovation with Quality Location”  In _X_ , 2022  URL: https://xmode.io/
  * [16] OpenStreetMap contributors  “Planet dump retrieved from https://planet.osm.org ”, https://www.openstreetmap.org, 2017
  * [17] Shiva Pourfalatoun and Erika E Miller  “Effects of covid-19 pandemic on use and perception of micro-mobility”  In _Available at SSRN 4113031_ , 2022
  * [18] Pitchaya Sakamanee, Santi Phithakkitnukoon, Zbigniew Smoreda and Carlo Ratti  “Methods for inferring route choice of commuting trip from mobile phone network data”  In _ISPRS International Journal of Geo-Information_ 9.5 MDPI, 2020, pp. 306
  * [19] Tejas Santanam, Anthony Trasatti, Pascal Van Hentenryck and Hanyu Zhang  “Public Transit for Special Events: Ridership Prediction and Train Optimization”  In _arXiv preprint arXiv:2106.05359_ , 2021
  * [20] Nilufer Sari Aslam, Tao Cheng and James Cheshire  “A high-precision heuristic model to detect home and work locations from smart card data”  In _Geo-spatial Information Science_ 22.1 Taylor & Francis, 2019, pp. 1–11
  * [21] Stefan Schönfelder et al.  “Analysis of commute Atlanta instrumented vehicle GPS data: Destination choice behavior and activity spaces”  In _Arbeitsberichte Verkehrs-und Raumplanung_ 350 ETH, Eidgenössische Technische Hochschule Zürich, IVT, Institut für …, 2006
  * [22] Anthony Trasatti  “Data-Driven Design of On-Demand Multimodal Transit Systems”, 2022
  * [23] Yang Xu et al.  “Understanding aggregate human mobility patterns using passive mobile phone location data: a home-based approach”  In _Transportation_ 42.4 Springer, 2015, pp. 625–646
  * [24] Xiping Yang et al.  “Understanding the spatial structure of urban commuting using mobile phone location data: a case study of Shenzhen, China”  In _Sustainability_ 10.5 MDPI, 2018, pp. 1435
  * [25] Wenbin Yao et al.  “Analysis of Key Commuting Routes Based on Spatiotemporal Trip Chain”  In _Journal of Advanced Transportation_ 2022 Hindawi, 2022
  * [26] Lei Zhang et al.  “Interactive covid-19 mobility impact and social distancing analysis platform”  In _Transportation Research Record_ SAGE Publications Sage CA: Los Angeles, CA, 2020, pp. 03611981211043813
  * [27] Xinhuan Zhang et al.  “Research on the method of travel area clustering of urban public transport based on Sage-Husa adaptive filter and improved DBSCAN algorithm”  In _PLoS one_ 16.12 Public Library of Science San Francisco, CA USA, 2021, pp. e0259472
  * [28] Yingjie Zhang, Beibei Li and Ramayya Krishnan  “Learning individual behavior using sensor data: The case of global positioning system traces and taxi drivers”  In _Information Systems Research_ 31.4 INFORMS, 2020, pp. 1301–1321