# Dedicated beam position monitor pair for model-independent lattice
characterization at NSLS-II

Yongjun Li<EMAIL_ADDRESS>Kiman Ha  Danny Padrazo  Bernard Kosciuk  Belkacem
Bacha  Michael Seegitz  Robert Rainer  Joseph Mead  Xi Yang  Yuke Tian  Robert
Todd  Victor Smaluk  Brookhaven National Laboratory, Upton, NY, 11973, USA
Weixing Cheng  Argonne National Laboratory, Lemont, IL, 60439, USA

###### Abstract

This paper reports recent lattice characterization results obtained at the
National Synchrotron Light Source II (NSLS-II) storage ring, conducted without
reliance on a lattice model. A pair of beam position monitors (BPMs) with
bunch-by-bunch (B$\times$B) resolution, were recently installed in a section
of the storage ring free of magnetic fields. The new BPM pair measured the
beam, or bunch’s transverse Poincaré map precisely after the beam was excited.
Linear one-turn-matrices (OTM) were then derived, and from these, the
4-dimensional coupled Twiss parameters were extracted at the locations of the
BPM pair. By normalizing beam oscillation amplitudes with the Twiss
parameters, the global action-variables were obtained. These action-variables
facilitated the measurement of the local Twiss parameters observed by other
BPMs independent on lattice model. This method is general, and particularly
useful in certain scenarios such as a round beam mode in a diffraction-limited
light source ring. We applied it to assess both weakly and strongly coupled
lattices at the NSLS-II ring. Through analysis of the strongly coupled
lattice, the quadrupole tilt errors were estimated to be less than 400
$\mu$rad. Utilizing the BPMs’ B$\times$B resolution, for the first time we
observed the variations of the linear lattice along a long bunch-train.

###### keywords:

Lattice characterization, model-independent method, coupling

##  1 Introduction

Model-independent lattice characterization methods using multi-turn data in
storage rings are well-documented across numerous studies, including ref. [1,
2, 3, 4]. This method is particularly useful when precise lattice models are
unavailable. For example, some extremely low-emittance light source rings will
operate in a round beam mode with on-resonance tunes to extend beam lifetime,
as noted by [5, 6]. In such cases, strong linear coupling can result from
unidentified random errors. Another scenario involves the obvious distortion
of the design lattice due to magnet imperfections and/or the influence of
insertion devices (IDs). Although ID fields can be taken into account at the
lattice design stage, during routine operations, their actual field can
dynamically vary due to changes in their gaps, phase adjustments and
compensation coils. Therefore, design models might not be able to accurately
represent a machine’s real-world operation. Given that at large user
facilities, their storage rings often host several tens of IDs, a model-
independent lattice characterization might be considered under these
conditions.

A direct method to characterize the linear lattice involves reconstructing the
one-turn-matrix (OTM) using multi-turn data captured by BPMs [7]. Assuming the
transfer matrix between a pair of the BPMs is known, with consecutive turn-by-
turn (T$\times$T) data, the Poincaré map can be constructed to fit OTMs.
However, this assumption is only approximately valid due to unidentified
random errors in the presence of real magnets, except when the BPM pair is
inside field-free region. In field-free regions the transverse momenta
(trajectory slopes) of the beam can be precisely determined using the formula:

$p_{x}\approx x^{\prime}=\frac{x_{2}-x_{1}}{L},$ (1)

where $x_{1,2}$ are the readings from the upstream and downstream BPMs at the
same turn, and $L$ is the distance between them. With well-calibrated and
aligned BPMs, the projected trajectories in the phase space can be accurately
reconstructed. To exploit this capability, a BPM pair was recently installed
in the cell 1 straight section of the NSLS-II ring. This section is reserved
for beam diagnostics and will remain free of magnets for the foreseeable
future.

In what follows, we summarize some recent lattice measurement results
utilizing the new BPM pair: Sect. 2 introduces the new BPMs’ electronics and
compares their T$\times$T resolution performance with other regular BPMs. In
Sect. 3, the lattice characterization results for both weak and strongly
coupled lattices are presented. Through fitting the measurement results, the
NSLS-II ring’s quadrupole tilt errors were estimated, for the first time,
based on a beam test. A preliminary observation of the non-linearity was also
presented. Sect. 4 reports the measurements of each individual bunch’s Twiss
parameters in a long bunch-train using their bunch-by-bunch (B$\times$B)
resolution. A brief summary is given in Sect. 5.

##  2 BPMs and their electronics

The recently installed BPM pair is situated within the cell 1 straight
section, separated by 5.14 $\mathrm{m}$. Their transverse cross-section
geometry is the same as our existing regular BPMs as detailed by [8]. Besides
regular functionalities, these BPMs feature advanced electronics capable of
providing B$\times$B resolution, allowing detailed exploration of the beam
dynamics at the bunch level. The electronics of these BPMs utilize a Radio
Frequency System-on-Chip (RFSoC) Field Programmable Gate Array (FPGA). The
RFSoC FPGA includes eight 14-bit Analog-to-Digital Converters (ADCs) with 5
Giga Samples Per Second (GSPS), directly sampling the BPMs’ four-button
signals. The BPM signals are synchronized with the 500 MHz RF reference clock
to align to each bunch. The system’s chassis incorporates a Xilinx ZCU208 FPGA
evaluation board and an RF breakout Balun board to interface with the single-
end analog input/output signals. Additionally, a custom-designed four-channel
RF switchboard and a hybrid Analog Front-End (AFE) module are employed for
beam signal conditioning. A linear AC/DC power supply minimizes amplifier
noise and an external Double Data Rate $4^{th}$-generation (DDR4) memory block
is integrated for fast, extensive data capture. Beam data is relayed to Input-
Output Controllers (IOCs) and client computers via the Experimental Physics
and Industrial Control System (EPICS). The functionalities of the electronics
firmware and signal processing are illustrated in Fig. 1. Further details on
the BPM electronics are discussed in ref. [9].

Figure 1: (Colored) Schematic diagrams of the B$\times$B BPM electronics. Top:
Diagram of the FPGA firmware functionalities. Bottom: Diagram of the FPGA
signal processing.

The performance of the new BPMs’ T$\times$T data acquisition was compared to
regular BPMs under identical conditions. T$\times$T data from a new BPM
(B$\times$B-2) and its nearest regular counterpart (BPM C01-1) are shown in
Fig. 2. The T$\times$T data prior to the $21^{st}$ turn, when the beam was not
kicked and was expected to remain stable. In this plot, the new BPM exhibited
significantly less fluctuation ($\approx 4\,$\mu\mathrm{m}$$) than its regular
counterpart ($\approx 22\,$\mu\mathrm{m}$$). In principle, T$\times$T
fluctuations can generally be smoothed out by averaging over multiple turns
for either the slow (10 Hz) or fast (10 kHz) orbit control. This average
results in a minimized impact on routine machine operations. The new
electronics, however, are particularly advantageous when precise T$\times$T
data is needed to measure the lattice functions. Additionally, the new BPM
pair’s B$\times$B resolution is approximately 2 to 3 times less precise than
the regular BPM’s T$\times$T resolution, and will be discussed in Sect. 4.

Figure 2: (Colored) Comparison of the T$\times$T readings of a new B$\times$B
BPM (B$\times$B-2) and its nearest regular counterpart, BPM (C01-1), when the
beam was excited at the $21^{st}$ turn. The new BPM’s noise fluctuation is
significantly smaller than the regular BPM before beam excitation.

##  3 Measurement of generally coupled linear lattices

###  3.1 Introduction to method

In a generally coupled linear lattice, beam Betatron oscillation is observable
by BPMs following excitation with the pulse kickers, referred to as pingers.
The turn-by-turn beam trajectory, as read by a BPM, can be mathematically
represented as:

$\begin{bmatrix}x_{i}\\\ y_{i}\end{bmatrix}=\begin{bmatrix}A_{1,x}\cos(i\cdot
2\pi\nu_{1}+\phi_{1,x})+A_{2,x}\cos(i\cdot 2\pi\nu_{2}+\phi_{2,x})+x_{co}\\\
A_{1,y}\cos(i\cdot 2\pi\nu_{1}+\phi_{1,y})+A_{2,y}\cos(i\cdot
2\pi\nu_{2}+\phi_{2,y})+y_{co}\end{bmatrix},$ (2)

where $\nu_{1,2}$ are the tunes of two oscillation modes, $i$ is the index of
turns, $\phi_{(1,2),{(x,y)}}$ are the initial phases, and $x_{co}$ and
$y_{co}$ represent the static closed orbits. The oscillation amplitudes are
given by:

$A_{(1,2),(x,y)}=\sqrt{2J_{1,2}\beta_{(1,2),(x,y)}(s)},$ (3)

where $J_{1,2}$ are the linear action-variables associated with the initial
excitation, which are global beam parameters consistent across all BPMs,
$\beta_{(1,2),(x,y)}(s)$ are the Twiss functions for modes $1$ or $2$ at the
location of $s$ in the horizontal $x$ or vertical $y$ planes, respectively. In
Eq. (2), the radiation damping and beam decoherence are ignored. A generalized
Courant-Snyder parameterization developed by Ripken et al., [10, 11, 12], is
adopted in Eq. (2). Alternative parameterization methods were available, such
as [13, 14, 15, 16], but were not applied in this study.

In our studies, we employed the harmonic analysis method [17] to determine the
Twiss $\beta$ functions. Utilizing the spectrum analysis techniques, such as
the Fast Fourier Transformation (FFT) or the Numerical Analysis of Fundamental
Frequencies (NAFF), the fractional tune $\nu_{1,2}$ can be obtained first. The
amplitudes of two harmonics (modes) in each plane were computed using their
respective cosine and sine components as follows:

$\begin{array}[]{ccc}C_{(1,2),(x,y)}&=&\sum_{i=1}^{N}(x,y)_{i}\cdot\cos(2\pi\nu_{1,2}i)\\\
S_{(1,2),(x,y)}&=&\sum_{i=1}^{N}(x,y)_{i}\cdot\sin(2\pi\nu_{1,2}i)\end{array}.$
(4)

By substituting Eqs. (2) into Eqs. (4), and averaging over a sufficient number
of turns (samples) $N$, the results are:

$\begin{array}[]{rrr}C_{(1,2),(x,y)}&\approx&\frac{N}{2}\sqrt{2J_{1,2}\beta_{(1,2),(x,y)}}\cos(\phi_{(1,2),(x,y)})\\\
S_{(1,2),(x,y)}&\approx&\frac{N}{2}\sqrt{2J_{1,2}\beta_{(1,2),(x,y)}}\sin(\phi_{(1,2)(x,y)})\end{array}.$
(5)

The amplitudes $A_{(1,2),(x,y)}$ of the oscillation observed at the BPM are
then calculated using:

$A_{(1,2),(x,y)}=\sqrt{2J_{1,2}\beta_{(1,2),(x,y)}}=\frac{2}{N}\sqrt{C_{(1,2),(x,y)}^{2}+S_{(1,2),(x,y)}^{2}},$
(6)

and the initial phases are determined by:

$\phi_{(1,2),(x,y)}=\arctan\left(\frac{S_{(1,2),(x,y)}}{C_{(1,2),(x,y)}}\right),$
(7)

where the sign of the phase $\phi$ depends on the signs of $C$ and $S$,
ensuring that the phase increases monotonically along the longitudinal
direction. This harmonic analysis approach (Eqs. (4)-(7)) can be conducted
using other methods such as the NAFF technique [18, 19].

The amplitudes $A$ in Eq. (6) are intertwined with the global action-variables
$J_{1,2}$, and the position-dependent Twiss functions
$\beta_{(1,2),(x,y)}(s)$. To determine the local $\beta$-functions, it is
essential to calibrate the action-variables, $J$, first. For weakly-coupled
linear lattices with an established model, such calibration can usually be
achieved by scaling the measured squared amplitudes, $A_{j}^{2}$ (with $j$ as
the index of BPMs) to match the design lattice model’s $\beta_{j}$ values.
This scaling approach, however, may introduce systematic errors in some
scenarios. One such example arises when strong coupling is introduced from
random unidentified errors with a near resonance tune. In this case, precise
lattice models are not always available. To address this challenge, it is
necessary to absolutely measure the $\beta$-functions at at least one BPM
location, and then extract the global action-variables with Eq. (6).

With the installation of the new BPM pair, the consecutive $4\times(N-1)$
phase space trajectories observed by the BPMs was used to derive the one-turn-
matrices $M$ as shown:

$\begin{pmatrix}x_{2},x_{3},\cdots,x_{n}\\\ p_{x,2},p_{x,3},\cdots,p_{x,n}\\\
y_{2},y_{3},\cdots,y_{n}\\\
p_{y,2},p_{y,3}\cdots,p_{y,n}\end{pmatrix}=M\begin{pmatrix}x_{1},x_{2},\cdots,x_{n-1}\\\
p_{x,1},p_{x,2},\cdots,p_{x,n-1}\\\ y_{1},y_{2},\cdots,y_{n-1}\\\
p_{y,1},p_{y,2}\cdots,p_{y,n-1}\end{pmatrix},$ (8)

where momenta $p_{x,y}$ were computed with Eq. (1) and $M$ is a $4\times 4$
matrix, with N representing the number of turns. From $M$, the Ripken Twiss
parameters at the BPM locations can be obtained, enabling the determination of
the global action variables $J$ via Eq. (6). These variables are then used to
normalize amplitudes measured by other BPMs. This methodology is generally
applicable to both weakly and strongly coupled linear lattices, therefore, no
detailed accelerator model is needed.

###  3.2 Measurement of a weakly coupled lattice

First we measured the $\beta$-functions and tunes for a weakly coupled
lattice, and compared against the design model. Given a weakly coupled
lattice, it is unnecessary to use this method, because the amplitudes can be
quite accurately normalized with the existing design model $\beta$-functions
as explained previously. However, we can take advantage of this to implement a
beam-based calibration for the new BPM pair, similar to the Linear Optics from
Closed Orbit (LOCO) technique [20].

A short bunch-train consisting of 10 consecutive buckets and carrying a charge
of 1.5 $\mathrm{m}\mathrm{A}$ was stored in the storage ring, which was
configured with a well-corrected linear lattice. Upon excitation, the 10
bunches began to perform coherent Betatron oscillations. With a T$\times$T
data set of 512 turns, the one-turn-matrix $M$ observed at the upstream BPM
was determined using Eq. (8):

${}M=\begin{pmatrix}[r]1.370\pm 0.002&4.608\pm 0.005&-0.020\pm 0.001&-0.221\pm
0.004\\\ -0.537\pm 0.001&-1.075\pm 0.002&-0.030\pm 0.001&-0.020\pm 0.001\\\
-0.008\pm 0.002&-0.171\pm 0.005&1.852\pm 0.003&5.751\pm 0.008\\\ -0.012\pm
0.001&0.025\pm 0.002&-0.843\pm 0.001&-2.079\pm 0.003\end{pmatrix}.$ (9)

The determinant of this matrix, $\det{M}=0.993\pm 0.001$, represents the
closest measurement to the symplectic condition ($\det{M}=1$) ever achieved at
the NSLS-II. From the matrix Eq. (9), a set of coupled Ripken Twiss functions
($\alpha,\beta,\gamma$) and two eigen-tunes (fractional parts) were extracted
and compared against their design values as shown in Tab. 1 and illustrated in
Fig. 3. The fractional tunes extracted from the T$\times$T data with the NAFF
technique closely align with those computed from the OTM. Although the
$\beta$-functions, dispersions and phase-advances were simultaneously
corrected prior to our measurement, each parameter retains has some residual
deviation. Therefore, the tune was not perfectly corrected to the design
value. Given the weak coupling in the lattice, the $\beta_{1x,2y}$ functions
degenerate into the commonly used, uncoupled $\beta_{x,y}$ functions.
Simultaneously, the minor values of $\beta_{1y,2x}$ confirm the coupling is
weak. Additionally, two new BPM roll errors were found negligible, and their
gains were calibrated as $g_{x}\approx g_{y}\approx 0.978$, as defined in the
LOCO algorithm [20]. The BPM calibration was based on comparisons with the
scaled action-variables averaged over 180 regular BPMs with a fitted model.
Note that those regular BPMs were previously calibrated with the LOCO
technique and T$\times$T data [21].

Table 1: Measured and design Ripken $\beta$-functions at the upstream BPM  name (unit) | mean | std. | design
---|---|---|---
$\beta_{1,x}$ ($\mathrm{m}$) | 4.659 |  $\pm$0.005 | 4.761
$\beta_{2,y}$ ($\mathrm{m}$) | 5.797 |  $\pm$0.008 | 5.765
$\beta_{1,y}$ ($\mathrm{m}$) | 0.007 |  $\pm$0.0003 | 0
$\beta_{2,x}$ ($\mathrm{m}$) | 0.004 |  $\pm$0.0001 | 0
$\alpha_{1,x}$ | 1.234 |  $\pm$0.0017 | 1.258
$\alpha_{2,y}$ | 1.984 |  $\pm$0.0031 | 1.981
$\alpha_{1,y}$ | 0.001 | $\pm 7\times 10^{-5}$ | 0
$\alpha_{2,x}$ | 0.003 | $\pm 1\times 10^{-4}$ | 0
$\nu_{1,OTM}$ | 0.2264 | $\pm 8\times 10^{-5}$ | 0.2213
$\nu_{2,OTM}$ | 0.2682 | $\pm 5\times 10^{-5}$ | 0.2598
$\nu_{1,NAFF}$ | 0.2264 | $\pm 8\times 10^{-5}$ | 0.2213
$\nu_{2,NAFF}$ | 0.2682 | $\pm 4\times 10^{-5}$ | 0.2598
Figure 3: (Colored) Comparison of the measured $\beta$-functions (dots with
error-bars) to the design model (dashed lines). The error bars, detailed in
Tab. 1, represent the means and standard deviations across 11 snapshots. Here,
$s=0$ $\mathrm{m}$ is the location of the injection point.

The phase space trajectories observed by the upstream BPM are illustrated in
Fig. 4. These trajectories are plotted alongside those predicted by the design
model comparative analysis. Notably, in the vertical plane, obvious smearing
was observed. Part of the smear might possibly be attributed to not-well-
corrected non-linearity of the lattice because of lack of independent
tunability of sextupoles powered in series [22]. In the presence of non-
linearity, to accurately extract the linear OTM, a sufficient number of data
points is needed. In our studies, we generally use 512 data points to mitigate
the nonlinear effect before Betatron oscillation is damped by the radiation
loss (the damping time: 20,800 turns [23]) and beam decoherence.

Figure 4: (Colored) Comparison of the measured 512 turns Poincaré maps (blue
dots) and the design model (oranges dots) as observed by the upstream BPM.

Minor deviations from the symplectic condition in $M$ could result from
factors such as the combination of BPM imperfections, radiation damping,
decoherence. A process known as “symplectification” was applied to the matrix.
In our case, the Hamilton-Cayley theorem [24] was applied to approximate $M$
to a linear Lie generator, $f_{2}=C_{abcd}x^{a}p_{x}^{b}y^{c}p_{y}^{d}$. Here
the indices $a,b,c,d$ are non-negative integers with a sum $a+b+c+d=2$, and
$C_{abcd}$ are the coefficients of monomials. Based on the Lie generator
$f_{2}$, another set of Twiss parameters and tunes were computed. The results,
however, were nearly identical to those derived directly with $M$.
Consequently, we concluded that the symplectification process could be
omitted.

The action-variables $J_{1,2}$ were determined using Eq. (6) as follows:

$J_{1}=0.0117\pm 1.8\times 10^{-5}\,$\mathrm{m}$,\;J_{2}=0.0104\pm 5.7\times
10^{-5}\,$\mathrm{m}$.$ (10)

Since the action-variable seen by all BPMs remains consistent, it was used to
normalize the amplitudes to determine the $\beta$-functions across the entire
ring as shown in Fig. 5. It is important to note that the measurement of the
Betatron phase, as demonstrated in the two bottom subplots, using Eq. (7),
neither needs a design model, nor depends on the BPM gains. However, a precise
alignment of all BPM data to the same turn and/or bunch is essential to
measure their phases correctly.

Figure 5: Comparison of the measured $\beta$-functions and phases (orange
lines with error-bars) against the design model (blue dotted-lines.), The
phases were normalized with $2\pi$ and then displayed with respect to a
modulus of 1.

###  3.3 Measurement of a strongly coupled lattice

In this section, we applied our method to a strongly coupled linear lattice in
which a lattice model is unavailable. Initially, all skew quadrupole
correctors were turned off and the tune was adjusted such that it was near the
difference resonance $\nu_{1}-\nu_{2}=17$. In this configuration, linear
coupling mainly arises from random and unidentified quadrupole tilts and
vertical orbit displacements through sextupoles, and no precise lattice model
is available. Since the beam orbit was already aligned to the centers of the
magnets using a beam-based alignment technique [25], the contribution from the
quadrupole tilts is considered dominant as explained in ref. [26]. Therefore,
by analyzing the quadrupole tilts with measured lattice functions, we can
estimate the overall misalignment of magnets and compare it with survey
results. As mentioned previously, the BPM gains and roll errors were pre-
calibrated and then excluded from the following coupling measurements.

From the T$\times$T data collected at 182 BPMs, a distinct beat pattern, due
to a strong linear coupling, can be observed as illustrated in Fig. 6. With
the absolute calibrated $\beta$-functions at the BPM pair, as explained in the
previous section, the four coupled $\beta$-functions seen by the upstream BPM
are: $\beta_{1x}=3.765\pm 0.017\,$\mathrm{m}$$, $\beta_{2y}=3.984\pm
0.035\,$\mathrm{m}$$, $\beta_{1y}=0.926\pm 0.016\,$\mathrm{m}$$,
$\beta_{2x}=0.948\pm 0.028\,$\mathrm{m}$$. The two global action-variables can
then be extracted as: $J_{1}=0.0113\pm 1.1\times
10^{-4}\,$\mathrm{m}$,\;J_{2}=0.0028\pm 1.4\times 10^{-4}\,$\mathrm{m}$$. The
coupled Twiss functions seen by the other 180 BPMs were determined by scaling
their amplitudes $A_{(1,2),(x,y)}$ with the actions $J_{1,2}$. The measured
lattice functions are illustrated in Fig 7.

Figure 6: An out-of-phase beat pattern as seen from the horizontal and
vertical T$\times$T data under strong linear coupling. Figure 7: (Colored)
Coupled Ripken $\beta$-functions observed by BPMs in a strongly coupled linear
lattice, with a near resonance tune $\nu_{1,2}=33.222/16.212$.

In a strongly coupled lattice, two linear action-variables projected onto two
planes can be calibrated independently, providing a method for self-
consistency checks. For example, $J_{1}$ can be calibrated either in the
horizontal or vertical plane as follows:

$J_{1}=\frac{A_{x}^{2}}{2\beta_{1,x}}=\frac{A_{y}^{2}}{2\beta_{1,y}},$ (11)

and similarly for $J_{2}$. The self-consistency check is illustrated in Fig.
8.

Figure 8: Self-consistency check of linear action-variables measured in two
planes, respectively. Each snapshot represents an independent measurement.

In Fig. 8, the action $J_{1}$ is approximately 3 times larger than $J_{2}$,
which results in a better signal-noise-ratio, thereby enhancing the self-
consistency measurements.

With the measured $\beta$-functions, we fitted a linear lattice model by
adjusting the tilt of all 300 quadrupoles. The model was constructed by using
an “approximate” Jacobian matrix:

$\textrm{J}_{\beta}=\left[\frac{\partial\beta_{(1,2),(x,y)}}{\partial\theta_{1}},\cdots,\frac{\partial\beta_{(1,2),(x,y)}}{\partial\theta_{n}}\right],$
(12)

here, $\beta$-functions are at the locations of BPMs, and
$\theta_{1,2,\cdots,n}$ are the $n$ quadrupole tilts. Lacking a precise linear
lattice model, this Jacobian matrix of dimension $720\times 300$ for
$180\times 4$ $\beta$-functions, and seen by 180 BPMs against 300 quadrupole
tilts, was computed with a “similar” near-resonance coupled lattice. This
similar lattice, which has the same tunes as the measured one, only represents
approximately the response of the $\beta$-function to the quadrupole tilt. The
near-resonance tune allows minor changes in quadrupole tilt to generate an
obvious coupled lattice change. The approximate Jacobian matrix was updated
iteratively through the fitting. After several iterations, a linear lattice
model which was reasonably consistent with the measurements, was obtained and
illustrated in Fig. 9. Convergence was generally achieved with the exception
of $\beta_{2x}$.

Some possible explanations for the imperfect convergence (particularly in
$\beta_{2x}$) may include: (1) some coupling contributions by the orbit’s
vertical offsets through sextupoles, indicating that merely fitting quadrupole
tilts might be insufficient; (2) the large number of fitting parameters
(quadrupole tilts in this case) relative to the objectives ($4\times$180
$\beta$-functions) may prevent optimal fitting of our model. Nevertheless,
based on the fitted model, we estimated the quadrupole tilt’s root-mean-square
(r.m.s.) to be around 400 $\mu$rad. As previously mentioned, part of the
coupling originates from the orbit offsets in sextupoles. Therefore, this
suggests that the actual quadrupole tilts at the NSLS-II ring should be better
than estimated. The specification of magnet tilts for installation is 200
$\mu$rad [23]. After the 10-years of operation and despite potential other
issues, such as settling of the floor and other infrastructure components, the
magnet alignment still maintains good accuracy. This is the first time we
estimated our quadrupole roll errors based on a beam test.

Figure 9: Fitting of a linear lattice model based on the measured
$\beta$-functions shown in Fig. 7.

Using the fitted model, the equilibrium eigen-emittances computed using the
SLIM technique [24, 26] were $\epsilon_{1,2}=1.6/0.5\,$\mathrm{n}\mathrm{m}$$.
This also allowed for estimation of local beam sizes throughout the entire
storage ring. Fig. 10 illustrates the transverse beam profiles at 15 short,
straight centers. The tilts of these profiles, represented by
$\tan(2\theta_{xy})=\frac{2\left<xy\right>}{\left<xx\right>-\left<yy\right>}$
were marked with dashed lines. Here $\left<xy\right>$, $\left<xx\right>$ and
$\left<xy\right>$ are the second moments of beam distributions. Variations in
beam size, and the presence of random tilts $\theta_{xy}$ result from the
randomness in quadrupole tilt distributions. The beam size information can be
used for computation of some critical beam parameters, such as beam lifetime
and undulator brightness.

Figure 10: Transverse beam profiles at 15 short straight centers (indexed at
the upright corners of subplots) by incorporating the fitted quadrupole tilts
into the lattice model. Dashed lines indicate the tilts of beam profiles
relative to the horizontal plane.

###  3.4 Comparison with a model-dependent measurement result

In this section, another model-dependent measurement of OTMs was implemented
across all 30 straight sections of the NSLS-II. Each section features a pair
of BPMs, flanking two sextupoles. In this configuration, trajectories in the
phase space can be reconstructed by incorporating the actual settings of the
magnets as illustrated in Fig. 11. By using the known magnet settings and the
readings of the two BPMs, the trajectories can be numerically resolved through
tracking simulations.

When nonlinear magnets such as sextupoles are positioned between two BPMs, and
only BPM readings are available, multiple trajectories can theoretically
exist. In practice, however, only the proximal trajectory is real and viable,
as others would extend beyond the dynamic aperture, and thus are not feasible.
Nevertheless, this method still relies on certain assumptions, including a
pre-established lattice model and known magnetic fields.

Using the same dataset as the previous section, the OTMs were reconstructed at
these BPM pairs. The symplectic quality of matrices, assessed by evaluating
their determinants, was found to be about 10% ($\left<\det{M}\right>=0.91$)
lower than those derived from new BPMs ($\det{M}=0.99$) as depicted in Fig.
12. The errors of $\beta$-functions and tunes derived from these non-
symplectic matrices can propagate and influence the characterization of the
full storage ring lattice.

Figure 11: Reconstruction of the proximal trajectory within an NSLS-II
straight section. The black dots represent the readings from a pair of BPMs
separated by two sextupoles. The black line depicts the reconstructed
trajectory in the horizontal plane, while the blue line represents the slope
$x^{\prime}$, which experiences jumps while passing through sextupoles. The
kicks received from the sextupoles are sub-$\mu$rad, therefore, the trajectory
appears almost linear between the BPM pair. Figure 12: Distribution of
$\det{M}$’s measured at NSLS-II straight sections where sextupoles intervene.
These $\det{M}$’s are not only significantly smaller than unity, but they also
exhibit large fluctuations. For context, the matrices derived from the new
BPMs show $\det{M}=0.9933\pm 0.0005$, illustrating better symplectic
integrity.

###  3.5 Observation of non-linearity

With the new BPM pair, the non-linearity of the lattice was observed more
accurately by projecting it into the normalized phase-space and action-angle
variables. Initially, the beam was excited with sufficiently small kicks to
keep its motion approximately linear, allowing for precise measurement of the
linear Twiss functions ($\beta$’s and $\alpha$’s) as shown in Tab. 1.
Subsequently, the beam was excited with gradually increasing kicks until non-
linearity became obvious. Notably, the new BPM pair allows for accurate
measurements of not only the $\beta$-functions, but $\alpha$-functions as
well. Through canonical transformations, we derived the $s$-independent
normalized Courant-Synder coordinates [27]
$\bar{x}=\frac{x}{\sqrt{\beta_{x}}},\,\bar{p}_{x}=\frac{\alpha_{x}x+\beta_{x}p_{x}}{\sqrt{\beta_{x}}}$,
and action-angle variables
$J_{x}=\frac{1}{2}\sqrt{\bar{x}^{2}+\bar{p}_{x}^{2}},\,\phi_{x}=\arctan\left(\frac{\bar{p}_{x}}{\bar{x}}\right)$
as illustrated in Fig. 13. The deformation of trajectories and the modulation
of action-variables with their conjugate angles served as measures of non-
linearity, which can potentially be used as objectives for online lattice
optimizations [28].

Figure 13: Visualization of lattice non-linearity in normalized horizontal
coordinates. Top: phase space trajectories in normalized Courant-Synder
coordinates. Bottom: Modulation of linear actions with their conjugate angle-
variables.

##  4 Bunch-by-bunch Twiss functions measurements

As mentioned previously in Sect. 2, this pair of BPMs can resolve individual
bunch Poincaré map with a 2 $\mathrm{n}\mathrm{s}$ bunch separation. This
section describes preliminary beam optics measurements that leverage this
capability. For this part of the investigation, 25 $\mathrm{m}\mathrm{A}$ beam
was stored in the storage ring, uniformly distributed in a bunch-train of
1200-buckets. As the NSLS-II storage ring has 1320-buckets, this leaves a gap
of 120 empty buckets. The beam was excited using pingers and a B$\times$B
dataset for all 1200-filled buckets was recorded over 1000 turns. By examining
their first turn data after pinger-excitation, we observed the profile of
initial excitation amplitudes along the bunch-train, which correlates well
with the pingers’ wave-forms as illustrated in Fig. 14. The pingers’ peaks
were aligned to reference bucket 1, i.e., the head of bunch-train. The flat
tops of the wave-forms were about 150 buckets (300 $\mathrm{n}\mathrm{s}$)
wide. For bunch trains exceeding this length, the T$\times$T data averaged
over the bunch train, post-excitation, is unsuitable for lattice measurements
due to incoherent bunch motions. Currently, the regular NSLS-II BPMs are
limited to the T$\times$T resolution; therefore, we rely on short bunch-
train’s T$\times$T data to characterize the linear lattice. It was noted that
the current B$\times$B resolution ($37$–$63\,$\mu\mathrm{m}$$, depending on
the bucket filling pattern and beam current) is approximately 2 to 3 times
less precise than the regular BPM’s T$\times$T resolution, as observed by the
fluctuations shown in Fig. 14. Despite this, the B$\times$B capability still
provides a potent diagnostic tool for investigating beam dynamics within bunch
trains, as explained below.

Figure 14: Initial excitation amplitudes along the bunch-train in the
horizontal and vertical planes. This profile correlates with the pinger’s
wave-forms.

Using the B$\times$B data, we reconstructed the Poincaré maps for each bunch
as shown in Fig. 15. Following reconstruction, their OTMs were fitted to
compute the Twiss parameter along the bunch train, depicted in Fig. 16.
Despite the relatively low stored beam current of 25 $\mathrm{m}\mathrm{A}$,
the linear Twiss functions have notable variation along the bunch-train,
particularly in the vertical plane. This variation suggests the existence of a
larger vertical than horizontal impedance. Considering the small vertical
aperture vacuum chambers used in dipole and ID sections, this result is
consistent with the existing NSLS-II infrastructure.

Should BPMs with the same B$\times$B capabilities be installed throughout the
entire storage ring, it would be feasible to re-construct the entire linear
lattice model for each individual bunch. Estimation of the quadrupolar wake
field distributions would also be possible. Additionally, another potential
application of the greater B$\times$B resolution would be visualization of
fast oscillating perturbations, such as the residual motion resulting from
imperfections in the injection transient.

Figure 15: Reconstructed phase space trajectories of the $30^{th}$ and
$250^{th}$ bunches at the upstream BPM. The $250^{th}$ bunch’s amplitudes
(measured with the greater B$\times$B capabilities) are slightly smaller than
the $30^{th}$ bunch’s due to its smaller initial excitation. From the
B$\times$B data, each individual bunch’s Twiss functions were extracted.
Figure 16: Variations in $\beta$-functions along a long bunch train, as seen
by the BPM pair. Here, the subscript “u” represents the upstream BPM and “d”
represents the downstream BPM. The vertical $\beta$-functions have a visible
variation, which might indicate larger vertical wake fields.

##  5 Summary

This paper summarizes recent lattice characterization efforts conducted at the
NSLS-II ring, facilitated by the installation of a new BPM pair. With the
installation location being free of magnetic fields, the BPM pair is
particularly useful when precise lattice models are unavailable, or when the
online machine significantly deviates from its design model. Through analysis
of a strongly coupled lattice, the quadrupole tilt errors were estimated.
Additionally, the unprecedented B$\times$B resolution of the BPM pair allowed
for observation of the linear lattice variation along a long bunch-train. It
is worth noting, however, with the greater instrumentation capabilities
provided by the BPMs discussed here, that more extensive studies are currently
taking place, but not covered in this manuscript. Some of these studies
include investigation into areas such as: collective beam effects with
different filling pattern configurations, as well as nonlinear dynamics
optimization through Poincaré maps analyses for large amplitude beam motions.

## Acknowledgements

We would like to extend our thanks to numerous NSLS-II colleagues for their
invaluable support, particularly Q. Shen, T. Shaftan, M. Capotosto, S.
Wilkins, S. Campbell, C. Yu and the accelerator coordination and operation
groups. Their contributions to the project management, hardware setup and beam
studies have been essential. This research was supported and funded by the
National Synchrotron Light Source-II’s Facility Improvement Projects (FIP),
and utilized resources of the NSLS-II, a U.S. Department of Energy (DOE)
Office of Science User Facility operated for the DOE Office of Science by
Brookhaven National Laboratory under Contract No. DE-SC0012704.

## References

  * [1] J. Irwin, C. Wang, Y. T. Yan, K. L. Bane, Y. Cai, F.-J. Decker, M. G. Minty, G. V. Stupakov, F. Zimmermann, Model-independent beam dynamics analysis, Physical review letters 82 (8) (1999) 1684.
  * [2] C.-X. Wang, Untangling mixed modes in model-independent analysis of beam dynamics in circular accelerators, Physical Review Special Topics-Accelerators and Beams 7 (11) (2004) 114001.
  * [3] X. Huang, J. Sebek, D. Martin, Lattice modeling and calibration with turn-by-turn orbit data, Physical Review Special Topics-Accelerators and Beams 13 (11) (2010) 114002.
  * [4] B. Riemann, P. Grete, T. Weis, Model-independent and fast determination of optical functions in storage rings via multiturn and closed-orbit data, Physical Review Special Topics-Accelerators and Beams 14 (6) (2011) 062802.
  * [5] M. Borland, L. Emery, R. Lindberg, V. Sajaev, Y. Sun, A. Xiao, Overview of lattice design and evaluation for the APS upgrade, ICFA Beam Dynamics Newsletter 71 (2017).
  * [6] C. Steier, A. Allézy, A. Anders, K. Baptiste, E. Buice, K. Chow, G. Cutler, R. Donahue, D. Filippetto, J. Harkins, et al., Status of the conceptual design of als-u, in: 9th International Particle Accelerator Conference, 2018, pp. 4134–4137.
  * [7] W. Fischer, Robust linear coupling correction with n-turn maps, Physical Review Special Topics-Accelerators and Beams 6 (6) (2003) 062801.
  * [8] O. Singh, B. Bacha, A. Blednykh, W. Cheng, L. Dalesio, C. Danneil, J. De Long, A. J. Della Penna Jr, K. Ha, Y. Hu, et al., NSLS-II BPM and fast orbit feedback system, IBIC13, Oxford, UK (2013).
  * [9] K. Ha, Y. Li, B. Bacha, D. Padrazo, V. Smaluk, Y. Tian, J. Mead, W. Cheng, NSLS-II bunch by bunch BPM development and beam operation, in: 15th International Particle Accelerator Conference Proc. IPAC’24, 2024, p. 2188.
  * [10] I. Borchardt, E. Karantzoulis, H. Mais, G. Ripken, Calculation of beam envelopes in storage rings and transport systems in the presence of transverse space charge effects and coupling, Zeitschrift für Physik C Particles and Fields 39 (3) (1988) 339–349.
  * [11] F. Willeke, G. Ripken, Methods of beam optics, in: AIP Conference Proceedings, Vol. 184, American Institute of Physics, 1989, pp. 758–819.
  * [12] V. Lebedev, S. Bogacz, Betatron motion with coupling of horizontal and vertical degrees of freedom, Journal of Instrumentation 5 (10) (2010) P10010.
  * [13] D. Edwards, L. Teng, Parametrization of linear coupled motion in periodic systems, IEEE Transactions on nuclear science 20 (3) (1973) 885–888.
  * [14] D. Sagan, D. Rubin, Linear analysis of coupled lattices, Physical Review Special Topics-Accelerators and Beams 2 (7) (1999) 074001.
  * [15] Y. Luo, Linear coupling parametrization in the action-angle frame, Physical Review Special Topics-Accelerators and Beams 7 (12) (2004) 124001.
  * [16] A. Wolski, Alternative approach to general coupled linear optics, Phys. Rev. ST Accel. Beams 9 (2006) 024001.  doi:10.1103/PhysRevSTAB.9.024001.
URL https://link.aps.org/doi/10.1103/PhysRevSTAB.9.024001

  * [17] J. Borer, C. Bovet, A. Burns, G. Morpurgo, Harmonic analysis of coherent bunch oscillations in LEP, Conf. Proc. C, 3rd European Particle Accelerator Conference (EPAC 92) 920324 (1992) 1082–1084.
  * [18] J. Laskar, C. Froeschlé, A. Celletti, The measure of chaos by the numerical analysis of the fundamental frequencies. application to the standard mapping, Physica D: Nonlinear Phenomena 56 (2-3) (1992) 253–269.
  * [19] P. Zisopoulos, Y. Papaphilippou, A. Streun, V. Ziemann, Beam Optics Measurements through Turn by turn Beam Position Data in the SLS, in: 4th International Particle Accelerator Conference, 2013, p. WEPEA067.
  * [20] J. Safranek, Experimental determination of storage ring optics using orbit response measurements, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 388 (1-2) (1997) 27–36.
  * [21] X. Huang, X. Yang, Linear optics and coupling correction with turn-by-turn BPM data, in: 6th International Particle Accelerator Conference Proc. IPAC’15, 2015, p. 698.
  * [22] Y. Li, D. Xu, V. Smaluk, R. Rainer, Nonlinear optics from hybrid dispersive orbits, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1060 (2024) 169032\.
  * [23] S. Dierker, NSLS-II preliminary design report, Brookhaven National Laboratory (2007).
URL https://www.osti.gov/biblio/1010602

  * [24] A. W. Chao, Lecture notes on topics in accelerator physics, Tech. Rep. SLAC-PUB-9574, SLAC National Accelerator Lab., Menlo Park, CA, USA (2002).
  * [25] G. Portmann, D. Robin, L. Schachinger, Automated beam based alignment of the als quadrupoles, in: Proceedings Particle Accelerator Conference, Vol. 4, IEEE, 1995, pp. 2693–2695.
  * [26] Y. Li, R. Rainer, Designing linear lattices for round beam in electron storage rings using the solution by linear matrices analysis, Phys. Rev. Accel. Beams 25 (2022) 040702.  doi:10.1103/PhysRevAccelBeams.25.040702.
URL https://link.aps.org/doi/10.1103/PhysRevAccelBeams.25.040702

  * [27] E. Courant, H. Snyder, Theory of the alternating-gradient synchrotron, Annals of physics 3 (1) (1958) 1–48.
  * [28] Y. Li, W. Cheng, L. Yu, R. Rainer, Genetic algorithm enhanced by machine learning in dynamic aperture optimization, Physical Review Accelerators and Beams 21 (5) (2018) 054601.


# gtfs2vec - Learning GTFS Embeddings for comparing Public Transport Offer in
Microregions

Piotr Gramacki<EMAIL_ADDRESS>0000-0002-4587-5586 Department of
Computational Intelligence, Wrocław University of Science and
TechnologyWrocławPoland ,  Szymon Woźniak<EMAIL_ADDRESS>0000-0002-2047-1649 Department of Computational Intelligence, Wrocław
University of Science and TechnologyWrocławPoland and  Piotr Szymański
<EMAIL_ADDRESS>0000-0002-7733-3239 Department of Computational
Intelligence, Wrocław University of Science and TechnologyWrocławPoland

(2021)

###### Abstract.

We selected 48 European cities and gathered their public transport timetables
in the GTFS format. We utilized Uber’s H3 spatial index to divide each city
into hexagonal micro-regions. Based on the timetables data we created certain
features describing the quantity and variety of public transport availability
in each region. Next, we trained an auto-associative deep neural network to
embed each of the regions. Having such prepared representations, we then used
a hierarchical clustering approach to identify similar regions. To do so, we
utilized an agglomerative clustering algorithm with a euclidean distance
between regions and Ward’s method to minimize in-cluster variance. Finally, we
analyzed the obtained clusters at different levels to identify some number of
clusters that qualitatively describe public transport availability. We showed
that our typology matches the characteristics of analyzed cities and allows
succesful searching for areas with similar public transport schedule
characteristics.

public transport timetable embeddings, typology of public transport offering,
unsupervised representation learning

††journalyear: 2021††copyright: acmlicensed††conference: 1st ACM SIGSPATIAL
International Workshop on Searching and Mining Large Collections of Geospatial
Data ; November 2, 2021; Beijing, China††booktitle: 1st ACM SIGSPATIAL
International Workshop on Searching and Mining Large Collections of Geospatial
Data (GeoSearch’21), November 2, 2021, Beijing, China††price: 15.00††doi:
10.1145/3486640.3491392††isbn: 978-1-4503-9123-8/21/11††ccs: Information
systems Geographic information systems††ccs: Computing methodologies Cluster
analysis††ccs: Computing methodologies Learning latent representations

Figure 1. Examples of regions with different types of public transport offer

##  1\. Introduction

Spatial data are available widely for many cities around the world. They are
used in solving urban-related tasks such as traffic flow prediction, house
price estimation, road network optimization, identification of accident
locations, managing transportation issues, and many others. They may also be a
valuable source of information in comparing so-called micro-regions in the
city. For this, however, it is necessary to have a method of representation of
those data. In this paper, we propose the first GTFS representation learning
approach which allows comparing regions in terms of the quality of public
transport offer. To confirm that the proposed representations are useful, we
perform in depth exploratory analysis of cluster obtained based on regions’
vector similarity. Using the obtained clusters, we formulate a public
transport offer typology and follow with examples of how the obtained vectors
can be used to search a database of micro-regions for regions with similar
public transport typology.

###  1.1. Representation learning in spatial data science

Term representation learning is used to describe a set of techniques for the
automated discovery of features from raw data. One of the methods from this
area is embedding. It can be defined as a task of finding an optimal function
$f:D\rightarrow\rm{I\\!R}^{d}$, which transforms data $D$ into $d$-dimensional
vector representation. The main advantage which is gained when translating
pure spatial data to latent space is the possibility to define a measure of
distance between objects. This makes comparing micro-regions in the city in a
certain spatial context possible. Finally, one can perform exploration in such
a new space, which may reveal previously unnoticed patterns or features
present in the city. Representation learning approach has been widely applied
to spatial data in recent years. Some of the works focus on proposing an
embedding method which is then tested on a variety of tasks (Du et al., 2018;
Xu et al., 2020; Wu et al., 2020), while other attempt to solve given tasks
and for this purpose, they apply representation learning techniques (Spruyt,
[n.d.]; Chang et al., 2018; Yao et al., 2018). They show an improvement when
compared to other methods used for solving such tasks (Chang et al., 2018; Du
et al., 2018; Yao et al., 2018; Xu et al., 2020; Wu et al., 2020). Moreover,
the obtained embedding spaces are often shown to be interpretable (Spruyt,
[n.d.]; Xu et al., 2020). Embedding methods for spatial data examined in this
review draw inspiration from other areas where the representation learning
approach is widely used. The majority of them is based on natural language
processing and word2vec (Mikolov et al., 2013) method (Chang et al., 2018; Du
et al., 2018; Yao et al., 2018). Others used existing methods for graph
embedding (Xu et al., 2020; Wu et al., 2020) or methods from the domain of
image processing (Spruyt, [n.d.]). This variety shows that there is a visible
and strong trend for applying representation learning for spatial data.

Public transport is picked up in the literature mostly in tasks related to
general delays estimation. Most recent methods apply embedding techniques to
improve the quality of delays estimating in public transport (Shoman et al.,
2020; Barnes et al., 2020). Delays were also used as a factor when comparing
regions around a city (Raghothama et al., 2016). The task of analyzing and
comparing regions in a city in terms of public transport availability is not
covered widely in the literature. Methods from the domain of network science
were successfully applied to this task (von Ferber et al., 2008), however,
this approach cannot be easily transferred to micro-regions comparison.
Utilization of travel times (Fayyaz S. et al., 2017) brings valuable insight
into the subject of public transport availability and can be transferred to
use micro-regions in place of stops. However, the computation costs of this
approach are significant, which opens an opportunity for proposing another
method that could be easier to apply to a large number of cities.

###  1.2. Public transport schedules

One of the urban types of data which to the best of the authors’ knowledge was
not deeply explored in the context of various machine learning or data mining
tasks is public transport schedules. Even though they are widely available in
a unified format for many cities around the world111https://transitfeeds.com/.
It is also obvious that public transport availability is an important aspect
when comparing regions in the city or comparing cities with each other. Hence
public transport quality is often one of deciding factors when choosing a
district to live in. Analysis of transport schedules in the sense of their
real suitability for city dwellers is not a trivial task. Therefore, providing
a method to learn meaningful vector representation of public transport
availability for a given micro-region seems to be important for the
development of systems with the aforementioned functionality. Ability to
analyze an entire city to find the areas which need an improvement in terms of
public transportation would be certainly beneficial for urban planners and
city authorities. This may also act as a way to compare current public
transport infrastructure between cities. This possibility may be used to
compare with cities considered as a ”gold standard” in the area of public
transport organization. Unfortunately, there is not much research conducted to
analyze an existing public transport infrastructure in terms of a static
analysis of public transport schedules. Many of them focus mainly on planning
optimal routes e.g. between different destinations. If they analyze existing
networks they often focus on the prediction of delays or calculation of the
estimated time of arrival (ETA). However, exploring similarities between
regions of the city or similarities between different cities is still a
relatively unexplored area.

General Transit Feed Specification (GTFS) is a unified format for static
public transport schedules introduced by Google (Google Inc., [n.d.]a). It was
first used to incorporate transit data to Google Maps and over the years
became de facto a standard format for transit feed data. It contains
information about stops, routes, and times or arrival for each stop on a
route. In addition to this basic set of information, the GTFS format allows
specifying the precise route which the vehicle takes. It can also store
additional information about a station such as levels, pathways, or fare
information.

###  1.3. Identification of micro-regions

A method of splitting a given area into micro-regions is a crucial aspect of
our work. There are many ways of achieving this goal. This section will go
through the most popular methods and will justify the choice made in this
paper.

The simplest solution is to divide a given area into regular shapes based on a
map projection. It is simple and intuitive but may result in some
inconveniences. Firstly, such division may depend on map scale, and therefore
may be inconsistent between different cities. Secondly, it heavily relies on
map projection which is inaccurate because the Earth is a sphere that cannot
be precisely mapped to a flat surface. Lastly, such division is not easily
repeatable as it may depend on the selected area’s borders.

Manual pointing of micro-region is the most accurate solution. It allows to
incorporate domain knowledge and creates a division suitable for a given task.
It makes it possible to divide a city based on different characteristics e.g.
on population density and/or incorporate administrative borders within the
city. An obvious disadvantage of such an approach is that it is labor-
intensive and specific to one single area. Therefore it is not suitable for
applications that will incorporate different regions of different cities, like
the one presented in this paper.

A spatial index makes it possible to divide a given area into grid cells,
identifiable by a unique index. Many spatial indexes work hierarchically,
meaning that each cell is composed of smaller cells. Examples of such systems
are Google’s S2 (Google Inc., [n.d.]b) based on rectangles and Uber’s H3 (Uber
Technologies Inc., [n.d.]) which uses hexagons. In this work, the H3 library
will be used to partition cities into micro-regions. The reasons behind this
decision are presented in the following section.

####  1.3.1. H3 hierarchical spatial index

Uber’s solution uses hexagons which is an effective way of space division. One
of its advantage over using rectangles or triangles is that all neighbors are
in the same distance. Another advantage is the hierarchical character of this
system. This allows selecting a resolution that is optimal for micro-region
identification and makes it possible to consider sub-regions for each micro-
region. Hierarchy in H3 is implemented as 16 levels of resolution. On level 0
there are 122 base cells covering an entire planet. The smallest hexagons on
level 15 have a diameter of length 1m.

For the task of regions embedding concerning public transport, two seem to be
a good fit - 8 and 9. They are respectively about 900m and 350m in diameter.
The bigger one covers an area roughly by main roads, whereas resolution 9 is
close to a single quarter of buildings. The difference is presented for a
small part of Wrocław in Figure 2.

Figure 2. H3 cells - comparison of resolution 8 and 9

Image presents a comparison between resolution 8 and 9 in Uber’s H3 spatial
index illustrated on a part of Wrocław. Each of bigger hexagon roughly
contains 7 smaller.

We analyzed how are stops distributed in regions on both resolutions. Using
resolution 8 means that there are regions with a lot of stops in them (over
30), whereas resolution 9 results in the majority of regions with no stop in
them. Based on that, we selected a resolution 8 further work.

The remainder of this paper is organized as follows: Section 2 presents a
proposed solution for features engineering and embedding model, Section 3
reports dataset utilized in this work and results of typology identification,
Section 4 summarizes whole work and presents future works. We release the
experimental code for this paper alongside obtained embeddings in a dedicated
Github repository222https://github.com/pwr-inf/gtfs2vec.

##  2\. Proposed solution

This section provides a detailed view of the proposed solution to the problem
of public transport accessibility representation and comparison between
cities. It covers a feature engineering process and normalization strategy.
After that, the architecture and learning process of an embedding model for
micro-regions representations learning is discussed in detail.

###  2.1. Features engineering

As presented in the first section, the GTFS format consists of multiple files
and is organized around a list of stop times for routes and stops. Such format
is impossible to be used directly in most ML/DL models. Therefore, we proposed
feature engineering work to prepare data in the correct format. This section
describes this process and justifies the applied solutions.

From a perspective of a single micro-region, two main public transport
accessibility metrics can be defined. Quantity and variety of public transport
available. Both were a part of the proposed solution.

####  2.1.1. Quantity related features

To describe the quantity of public transport available for a given micro-
region, we calculated the sum of trips in a given period. To perform those
calculations, we used the gtfs-kit library (Raichev, 2019). It provides
methods for the calculation of trips counts for each stop at a given time-
frequency. Those numbers were later summed up across an entire micro-region.

Sometimes multiple stops, which form the same route, are aggregated into a
single micro-region. As a result, the same trip will be counted multiple
times. However, this is an intentional assumption, because those features are
not meant to have a meaning of counting vehicles leaving a micro-region.
Instead, they describe in how many ways one can leave an area. The intuition
behind this is that having more options to get on a bus/tram inside a region
is better than just a single stop in the middle of an area.

####  2.1.2. Variety related features

Another feature which can describe public transport availability is the
variety of directions in which one can leave a region. It can be considered on
different levels of precision starting from world directions and ending on all
regions accessible via public transport with a given number of transfers. For
this thesis, an in-between solution was selected. It is described below.

The part of a route description in GTFS format is a trip_headsign column. In
practice, this is the destination of a route displayed on a vehicle. They are
in most cases either transportation hubs or line ends located most often on
the outer ring of the city. This is in a balance between the generality of
world directions and the specificity of each region connected via public
transport. Hence, to create a feature for a region representation, unique
head-signs were counted in a given period, resulting in the total number of
unique directions available from a region.

####  2.1.3. Time resolution for features and filtering

Both of the features mentioned above can be calculated for different time
windows resulting in a different number of features. For this work, one hour
window was selected resulting in 24 quantity and 24 variety features. This
resolution seems to be a reasonable choice because it presents changes
throughout the day and allows to interpret periods for example as rush hours
or nighttime.

After an analysis of timetables for night lines, we decided to remove them
from features of a region. Nighttime transportation availability is
significantly different from daytime, as the routes are often longer to
minimize the number of vehicles and drivers needed. Moreover, some additional
hubs were identified which are not active during the day. Therefore, only 6-22
hours were included as features for each region. As a result, each micro-
region was described with a vector $\mathbf{x}$ of length 34.

###  2.2. Normalization

The main purpose of our work was to propose a method of learning
representation of public transport accessibility described by features
mentioned in the previous section. As an embedding model, we used a neural
network, which will be described in the following section. Here we must note
that the above-mentioned features which are the main data in our task should
be normalized to be able to serve as an input to the neural network.

The values of features describing the number of trips are approximately 10
times larger. This may result in a neural network focusing mostly on those
features. Another observation is that values vary for different cities. This
was also taken into consideration when proposing a method of normalization.

To achieve the best performance of the embedding model, features will be
scaled to the $[0,1]$ range using min-max normalization technique, which is
expressed with the formula

(1)
$\mathbf{x}^{\prime}=\frac{\mathbf{x}-\min(\mathbf{x})}{\max(\mathbf{x})-\min(\mathbf{x})}.$

As values of each type (quantity and variety) should maintain relative values
during the whole day, normalization will be performed twice: once for all
columns with a sum of trips and once for all columns with destinations counts.

Normalization was conducted for all cities combined. This will ensure that
differences in public transport available between cities are maintained. This
approach should allow comparing the quality of public transport between
cities.

###  2.3. Embedding method

Autoencoder is a special type of neural network capable of learning lower- or
higher-dimensional embeddings of data. It is trained in an unsupervised
manner, which is crucial for this thesis.

Autoencoder network consists of two networks - encoder and decoder. The
encoder network can be treated as a function
$\phi:\mathcal{X}\rightarrow\mathcal{Z}$, and decoder as
$\psi:\mathcal{Z}\rightarrow\mathcal{X}$. Having input vector
$\mathbf{x}\in\mathcal{X}$, encoder network produces an embedding
$\phi(\mathbf{x})=\mathbf{z}$, which is then fed into decoder network to
obtain reconstructed input $\psi(\mathbf{z})=\mathbf{x^{\prime}}$. Both
networks are trained together using backpropagation algorithm, with the loss
function for a batch of inputs $\mathbf{X}$ defined as:

(2)
$\mathcal{L}=\frac{1}{N}\sum_{n=1}^{N}(\mathbf{X}_{n}-\mathbf{X^{\prime}}_{n})^{2},$

which is a mean squared error loss. When the training is complete, only the
encoder is used for inference.

Historically first applications of autoencoders focused on compressing an
input and discovering nonlinear correlations in data can be found in the work
(Kramer, 1991). It is possible to use nonlinear activation functions, and
smaller dimensionality of space $\mathcal{Z}$. In this work, we used the
opposite approach. We defined a sparser, 64-dimensional representation with an
architecture presented below.

  * •

Encoder:

    * –

$34\rightarrow 48$ neurons fully connected layer,

    * –

ReLU activation function,

    * –

$48\rightarrow 64$ neurons fully connected layer;

  * •

Decoder:

    * –

$64\rightarrow 48$ neurons fully connected layer,

    * –

ReLU activation function,

    * –

$48\rightarrow 34$ neurons fully connected layer;

This particular architecture was selected somewhat arbitrarily because in a
representation learning domain the goal is to extract representations from raw
data and only at the stage of exploratory analysis it is possible to determine
whether they are good or not. However, several intuitions are backing up this
choice of architecture. To begin with, the amount of data, which is available
for training is not big. The reason for that is the limited availability of
timetables data. A city is divided into approximately 300-500 micro-regions,
which means that there is not enough data to train a very deep model.
Therefore, we decided not to add any more layers to our model. Moreover, a
sparser representation was used to provide more dimensions to differentiate
types of public transport. Having considered that, a 64-dimensional embedding
layer and one hidden layer seem like a reasonable choice.

##  3\. Similarities detection

A final goal of analysis was to define a typology of public transport offer
types across multiple cities. To achieve that, we ran clustering multiple
times increasing the number of clusters at each step. At each step, we
analyzed obtained clusters to determine which number of clusters forms an
interpretable typology. The last step will be to describe an obtained
typology.

###  3.1. Datasets

In this research we used Open Mobility Data (MobilityData IO, [n.d.]) service
to gather GTFS files for multiple cities in Europe.333with an exception of
Zurich, which was downloaded from https://data.stadt-
zuerich.ch/dataset/vbz_fahrplandaten_gtfs To limit the number of cities we
defined the following inclusion criteria:

  * •

city population of at least 200,000 people,

  * •

at least 20 routes of public transport,

  * •

feed must be complete, meaning no public transport types should be missing,

  * •

GTFS file should not be for an entire country (those include intercity trains,
which may affect feature engineering process).

While processing feeds we encountered some issues, which resulted in an
additional exclusion criterion - a trip_headsign column in the trips table
must be correct, since it is used to calculate a number of directions
available from a region. As a result, we collected data for 48 cities, which
are listed in Table 1.

Table 1. List of cities included in research Aachen | Antwerp | Athens | Barcelona | Belgrade
---|---|---|---|---
Berlin | Bilbao | Brussels | Budapest | Bydgoszcz
Cologne | Dublin | Erfurt | Espoo | Florence
Genua | Gent | Hamburg | Helsinki | Karlsruhe
Krakow | Lille | Lublin | Madrid | Mannheim
Milan | Nantes | Nice | Oslo | Oulu
Palermo | Poznan | Prague | Radom | Szczecin
Tampere | Toulouse | Tricity (Gdansk, Gdynia, Sopot)
Turin | Valencia | Venice | Vienna | Vilnius
Warsaw | Wroclaw | Zagreb | Zurich | Leipzig

###  3.2. Agglomerative clustering

Having an embedding model trained, we calculated embeddings for all micro-
regions. The next step was to analyze the results and finally try to extract
valuable information from comparing the cities. This step will use
unsupervised methods like distance calculations in embedding space later used
in clustering task. As the task solved in this thesis is purely unsupervised
in nature, there is no other way to evaluate results other than analyzing
obtained space.

The clustering method requires a distance metric to be defined in obtained
embedding space. As a distance metric, cosine and euclidean metrics were
considered, following a common practice from literature. A cosine metric finds
similarities in direction from the origin, which would discard quantity
differences between values and focus on a general trend. This would work great
in a solution, which aims to find similarities between regions in terms of
public transport usage type because it discards differences in for example
number of trips. However, in this work, that those differences are important.
Therefore, a euclidean distance is used, since it extracts actual distances
between samples. It is calculated as

(3) $d(\mathbf{x},\mathbf{y})=||\mathbf{x}-\mathbf{y}||_{2}$

The last element of the proposed solution for exploratory analysis of obtained
embedding space is the clustering method based on a defined distance metric.
One of the main goals of this analysis was to propose a multi-level
classification in terms of public transport accessibility and quality. This
was a reason for choosing a hierarchical approach to clustering, namely an
agglomerative clustering.

This is an iterative method, which starts by placing all samples in separate
clusters. Then, at each step, two nearest clusters are merged, until a given
number of clusters is obtained or another stop criterion is met. To select
clusters to merge, Ward’s method (Murtagh and Legendre, 2014) was used,
because it minimizes intra-cluster variance. This feature is beneficial when
identifying a well-differentiated typology.

To better understand a character of obtained clusters we used two aggregated
features: a) sum_trips which represents the sum of trips from an entire day
(within 6-22 limits), b) directrions_whole_day which represents all directions
available during a day (from 6 to 22). They represent public transport
availability in a region from whole day perspective. This is a good measure to
justify differences between clusters from a distance.

On the first level of typology, we identified three major types of areas. They
were called suburban, mid-city, and hubs for this work. Clearly, those
clusters describe a general availability of public transport, where suburban
has relatively low availability and mid-city has high. Hubs are typical
outliers with very rich public transport offer.

Figure 3. First level typology visualized for Berlin. Figure 4. Second level
typology visualized for Berlin.

Figure 5. Line cities organized alongside a train line - Zagreb and Tricity
(Gdansk, Gdynia, Sopot)

###  3.3. Typology definition

The second level utilized the division of each type into three more. Each of
the level one types was divided differently. From suburban areas, we extracted
regions with high diversity and high quantity of available public transport
(respectively clusters 1 and 7). They can be treated as local transfer points
in an area with poorer public transport availability and well-connected parts
of suburbs (for example just one bus line but with a high frequency of trips).
The remaining regions from suburban-type have significantly worse public
transport availability.

The division in mid-city type was done similarly. Regions with a high
frequency of public transport but with limited variety were separated (cluster
3). Those are regions with very good availability of public transport. The
second type separated regions with a relatively big variety of public
transport (cluster 5). Those can serve as transportation hubs and main
transfer points in cities with worse quality of public transport network or as
local hubs in districts. The remaining cluster includes regions with city hubs
which separated differently. One of the clusters shows significantly better
public transport availability (cluster 2). It represents hubs with the highest
number of trips and directions. Two remaining clusters maintain only one of
them on a high level - either number of trips (cluster 6) or a number of
directions (cluster 8). The first one would mean regions with a dense metro
network. The other contains quite often train stations. Summarized typology
can be defined as:

  * •

suburban areas

    * –

cluster 4 \- regions with the least amount of public transport available

    * –

cluster 1 \- parts of suburbs with better variety of public transport - some
sort of entry points for suburbs

    * –

cluster 7 \- regions in suburbs with more trips but with small variety of
public transport

  * •

mid-city

    * –

cluster 0 \- mid-city regions with regular public transport availability

    * –

cluster 3 \- mid-city regions with a lot of trips every hour on limited number
of directions

    * –

cluster 5 \- regions with very good public transport availability, in cities
with worse public transport network they may serve as hubs, in other just as
regional hubs

  * •

hubs

    * –

cluster 2 \- very high intensity hubs

    * –

cluster 6 \- regions with a lot of trips, but limited directions, most often
regions with metro network in a city center

    * –

cluster 8 \- railway-type hubs with a lot of diversity and less trips per hour

Both levels of the typology were visualized on a map which is available
online444
https://kepler.gl/demo/map?mapUrl=https://dl.dropboxusercontent.com/s/qge8lyeinpmm7ud/keplergl_53g73fa.json.
The map includes three layers for each of the evaluated cities:

  * •

first level typology - with 3 clusters, as depicted for Berlin in Figure3;

  * •

second level typology - with 9 clusters, as depicted for Berlin in Figure4;

  * •

layer with numeric features from the input table, with number of trips or
destinations.

Figure 6. Imbalance in public transport availability between left- and right-
bank part of Warsaw

###  3.4. Similarities between cities

Our method allows for identification of regions with similar role in different
cities. We used both levels of typology to identify those similarities, where
the first level allows us to compare cities from a bird’s eye view. We used it
to compare the types of public transport in suburban areas in different
cities. First example is a so called ”line city” organized around a train
line. Two examples of Zagreb and Gdansk are presented in Figure 5. They both
have regions from mid-city cluster along a train line, and one or more hubs
close to main train stations. Parts of a city, which are further from a train
line have poorer public transport availability.

Our vector representations also allow finding imbalances in public transport
availability in different parts of a city. Figure 6 shows a view of Warsaw,
with a visible separation for left- and right-bank parts. Even though, both
parts are similar in size, a right-bank part of a city has worse public
transport availability. Same observations can be made for Vienna and Budapest
(which is visible on maps, which we provide a link to in this publication).

When analyzing a presence of regions assigned to the hubs cluster on first
level of typology, we identified three types of cities. First is the one
without those regions - smaller cities or the ones where public transport
stops are not do dense to form a regions of this characteristic (an example is
Poznan or Cologne). Second type have those hubs in city center, most often in
the neighborhood of a train station (an example can be Gdansk, Gdynia and
Zagreb in Figure 5). Final type is a city with multiple hub-like regions which
form bigger groups (see Berlin in Figure 3 or Warsaw in Figure 6 as well as
Madrid, Athens or Budapest which can be seen on a linked map). It represents
cities with the highest availability of public transport or with a very dense
public transport network in a city center.

Using the second level of an obtained typology, we were able to extract
regions with similar role across the cities. In this part, we present those
regions to show, that our method is efficient of matching regions extraction
between cities. We focused on regions in a city center, which can serve
different purposes in cities public transport network.

Figure 7. Regions with train stations inside are detected just based on public
transport availability. From left to right: Aachen, Berlin, Gent, Hamburg,
Karlsruhe and Tampere.

First of the identified region types is a region with railway station in it.
It was marked in this work as cluster 8 and in Figure 7 we present examples
from Aachen, Berlin, Gent, Hamburg, Karlsruhe and Tampere, where in each of
those cities this cluster is assigned to the region with a train station
inside. High availability of public transport in those regions is a
combination of a railway connections and public transport stops in a close
proximity.

Figure 8. Regions from a center of cluster with high public transport
availability and variety: Cracow, Berlin, Prague and Oslo.

Second type of region in a hub category is represented by cluster 2 and
combines regions with very high availability of public transport. In Figure 8
we presented 4 out of 5 regions from the center of cluster 2. They land in
Cracow, Berlin, Prague and Oslo and all contain regions with high number of
stops of public transport available and serve the role of transportation hubs.

Another type of region is a local hub associated with cluster 5. Our analysis
revealed, that those regions, when located further from a city center, serve a
role of local transportation hubs/transfer points. Examples are shown in
Figures 4 (Berlin) and 10 (Wroclaw). In Berlin they mostly match with U-Bahn
or railway stations and in Wroclaw they cover major junctions where multiple
public transport lines cross.

Figure 9. 2nd level of typology visualized for Wroclaw, regions in cluster 5
with local hubs are clearly discernible

The last example shows an analysis of a standing out city. Public transport
network in Madrid is very dense and therefore it stands out collections most
of the regions associated in cluster 6. Almost whole city is classified to a
type with high public transport availability. This is probably due to a very
dense subway network in this city. This shows, that our proposed method is
capable of not only identifying similar regions in different cities, but also
finding the cities which stand out in terms of public transport availability.

Figure 10. Madrid’s dense public transport network in the city center - an
(standing out) example of cluster 6 \- regions with a lot of trips, but
limited directions.

##  4\. Conclusions and future work

We proposed a solution to obtain vector representation of public transport
offer in GTFS format. Using exploratory cluster analysis we show that it
captures transport characteristics and allows identifying similarities between
regions. Thanks to using a wide selection of cities allowed for a good
quality, our multi-level typology definition can be used for comparisons
between cities and within them. Moreover, several unique functions of regions
were identified in terms of their role in the public transport network.
Finally we succesfully used obtained embeddings to search for regions with
similar public transport offers across evaluated cities, which a novel
contribution to the field.

In the future investigating other forms of normalization can be an interesting
direction. Usage of the local normalization approach, meaning that
normalization is conducted separately for each city. This approach should
allow for independence from the scale of available transportation in a city
and allow for the identification of regions with similar functions despite
different public transportation availability. At the cost of being able to
compare cities among themselves, this approach should allow comparing the
styles of public transport organizations in different cities.

## References

  * (1)
  * Barnes et al. (2020) Richard Barnes, Senaka Buthpitiya, James Cook, Alex Fabrikant, Andrew Tomkins, and Fangzhou Xu. 2020\.  BusTr: Predicting Bus Travel Times from Real-Time Traffic.  _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ (jul 2020), 3243–3251.  https://doi.org/10.1145/3394486.3403376 arXiv:2007.00882
  * Chang et al. (2018) Buru Chang, Yonggyu Park, Donghyeon Park, Seongsoon Kim, and Jaewoo Kang. 2018.  Content-Aware Hierarchical Point-of-Interest Embedding Model for Successive POI Recommendation. In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence_ , Vol. 2018-July. International Joint Conferences on Artificial Intelligence Organization, California, 3301–3307.  https://doi.org/10.24963/ijcai.2018/458
  * Du et al. (2018) Jiahong Du, Yujun Chen, Yue Wang, and Juhua Pu. 2018\.  Zone2Vec: Distributed Representation Learning of Urban Zones. In _2018 24th International Conference on Pattern Recognition (ICPR)_ , Vol. 2018-Augus. IEEE, 880–885.  https://doi.org/10.1109/ICPR.2018.8545376
  * Fayyaz S. et al. (2017) S. Kiavash Fayyaz S., Xiaoyue Cathy Liu, and Guohui Zhang. 2017\.  An efficient General Transit Feed Specification (GTFS) enabled algorithm for dynamic transit accessibility analysis.  _PLOS ONE_ 12, 10 (oct 2017), e0185333.  https://doi.org/10.1371/journal.pone.0185333
  * Google Inc. ([n.d.]a) Google Inc. [n.d.]a.  Reference — Static Transit — Google Developers.
https://developers.google.com/transit/gtfs/reference.  Accessed: 2021-04-05.

  * Google Inc. ([n.d.]b) Google Inc. [n.d.]b.  S2 Geometry.
https://s2geometry.io/.  Accessed: 2021-04-23.

  * Kramer (1991) Mark A. Kramer. 1991\.  Nonlinear principal component analysis using autoassociative neural networks.  _AIChE Journal_ 37, 2 (feb 1991), 233–243.  https://doi.org/10.1002/aic.690370209
  * Mikolov et al. (2013) Tomas Mikolov, Kai Chen, Greg S. Corrado, and Jeffrey Dean. 2013\.  Efficient Estimation of Word Representations in Vector Space.  http://arxiv.org/abs/1301.3781
  * MobilityData IO ([n.d.]) MobilityData IO. [n.d.].  Open Mobility Data.  https://openmobilitydata.org/.
Accessed: 2021-04-16.

  * Murtagh and Legendre (2014) Fionn Murtagh and Pierre Legendre. 2014.  Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?  _Journal of Classification_ 31, 3 (oct 2014), 274–295.  https://doi.org/10.1007/s00357-014-9161-z
  * Raghothama et al. (2016) Jayanth Raghothama, Vinutha Magal Shreenath, and Sebastiaan Meijer. 2016.  Analytics on public transport delays with spatial big data. In _Proceedings of the 5th ACM SIGSPATIAL International Workshop on Analytics for Big Geospatial Data, BigSpatial 2016_. Association for Computing Machinery, Inc, New York, New York, USA, 28–33.  https://doi.org/10.1145/3006386.3006387
  * Raichev (2019) Alex Raichev. 2019\.  GTFS Kit.  Github repository: https://github.com/mrcagney/gtfs_kit.
  * Shoman et al. (2020) Maged Shoman, Armstrong Aboah, and Yaw Adu-Gyamfi. 2020\.  Deep Learning Framework for Predicting Bus Delays on Multiple Routes Using Heterogenous Datasets.  _Journal of Big Data Analytics in Transportation_ 2, 3 (dec 2020), 275–290.  https://doi.org/10.1007/s42421-020-00031-y
  * Spruyt ([n.d.]) Vincent Spruyt. [n.d.].  Loc2Vec: Learning location embeddings with triplet-loss networks.  https://www.sentiance.com/2018/05/03/venue-mapping/.  Accessed: 2021-03-21.
  * Uber Technologies Inc. ([n.d.]) Uber Technologies Inc. [n.d.].  H3: Uber’s Hexagonal Hierarchical Spatial Index.
https://eng.uber.com/h3/.  Accessed: 2021-04-20.

  * von Ferber et al. (2008) C. von Ferber, T. Holovatch, Yu. Holovatch, and V. Palchykov. 2008.  Public transport networks: empirical analysis and modeling.  _European Physical Journal B_ 68, 2 (mar 2008), 261–275.  https://doi.org/10.1140/epjb/e2009-00090-x arXiv:0803.3514
  * Wu et al. (2020) Ning Wu, Xin Wayne Zhao, Jingyuan Wang, and Dayan Pan. 2020\.  Learning Effective Road Network Representation with Hierarchical Graph Neural Networks. In _Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. Association for Computing Machinery, 6–14.  https://doi.org/10.1145/3394486.3403043
  * Xu et al. (2020) Shuai Xu, Jiuxin Cao, Phil Legg, Bo Liu, and Shancang Li. 2020.  Venue2Vec: An Efficient Embedding Model for Fine-Grained User Location Prediction in Geo-Social Networks.  _IEEE Systems Journal_ 14, 2 (jun 2020), 1740–1751.  https://doi.org/10.1109/JSYST.2019.2913080
  * Yao et al. (2018) Zijun Yao, Yanjie Fu, Bin Liu, Wangsu Hu, and Hui Xiong. 2018.  Representing Urban Functions through Zone Embedding with Human Mobility Patterns. In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence_ , Vol. 2018-July. International Joint Conferences on Artificial Intelligence Organization, California, 3919–3925.  https://doi.org/10.24963/ijcai.2018/545