080003 2016 L. A. Pugnaloni C. M. Carlevaro, Instituto de Física de Líquidos y
Sistemas Biológicos,
La Plata, Argentina. 080003

We describe a simple experiment involving spheres rolling down an inclined
plane towards a bottleneck and through a gap. Results of the experiment
indicate that flow rate can be increased by placing an obstruction at optimal
positions near the bottleneck. We use the experiment to develop a computer
simulation using the PhysX physics engine. Simulations confirm the
experimental results and we state several considerations necessary to obtain a
model that agrees well with experiment. We demonstrate that the model exhibits
clogging, intermittent and continuous flow, and that it can be used as a tool
for further investigations in granular flow.

# Increasing granular flow rate with obstructions

Alan Murray  [inst1] Fernando Alonso-Marroquin[inst2] E-mail:
<EMAIL_ADDRESS>
(22 November 2015; 19 February 2016)

††volume: 8

99 inst1 SAE Creative Media Institute, Sydney, Australia. inst2 School of
Civil Engineering, The University of Sydney, NSW 2006, Australia.

##  1 Introduction

When does an obstruction placed near a bottleneck increase the flow rate of
discrete objects moving through the bottleneck? Answering this question is of
great utility on many scales. For example, the safe evacuation of pedestrians
moving through confined environments (train stations, stadiums, concert halls,
etc) in emergency situations can be of vital importance. And grain falling
through a silo, where efficient flow is necessary for production and clogging
is undesirable, is another example.

As stated by Magalhaes et al. in Ref. [1], ‘Granular materials are ubiquitous
either in nature or in industrial processes’ and so a fundamental
understanding of their motions is of intrinsic interest, both from physics and
engineering perspectives. The quantitative study of the flow of granular
materials has been performed for many decades [2], yet there is enormous scope
for many interesting and creative techniques to be developed: experimental,
theoretical and computational. Granular flow around an obstruction placed near
a bottleneck is of particular interest as it has many potential applications
and this has been investigated by many authors [1, 3, 5, 6, 7, 4].

Increases in flow rate achieved through placement of an obstruction near a
bottleneck have been reported by several investigators [8, 9, 10]. Another
commonly reported phenomenon is that of ‘clogging’, also described as the
formation of arches [4]. There are also many novel investigations of granular
flow. For example, Wilson et al. [3] consider granular flow of particles that
are completely submerged under water and Lumay et al. [7] investigate the flow
of charged particles that repel each other. Clearly, granular flow is a very
wide and active area of research.

Several authors have investigated the use of physics engines used in games, an
example of which can be found in the work by Carlevaro and Pugnaloni [15]. In
this paper we present a real time, 3D computer simulation using the PhysX
physics engine. The model can be used in the study of granular flow as it
exhibits continuous flow, intermittent flow, clogging, and has several
parameters that an investigator can control as the simulation is occurring.
The computer model is based on an experiment that uses spheres for the
particles and three different obstruction shapes. We use experimental data to
calibrate and validate our model. The model can then be used to determine how
the highest flow rates may be achieved both with regard to obstruction
position and to find optimal obstruction shapes.

The paper is organized as follows: we present the model in section 2 Then we
describe the experiments in two parts. Section III deals with flow rate
without obstructions, and section IV investigates flow rate with obstructions.
We discuss the limitations in section V, and conclude in section 6

##  2 The Model

We use PhysX [12], a freely available real time, 3D physics engine for
computations. PhysX is widely used in the games industry as, when necessary,
it favors speed and stability of computation over accuracy. Integration of the
equations of motion is done using an unconditionally stable, semi implicit
Euler scheme yielding algebraic equations that are solved using a progressive
Gauss-Seidel algorithm, whilst enforcing the Signorini conditions [16]. The
time step can be adjusted for greater accuracy at the cost of speed and hence,
real time performance of the engine. For this study, we left the time step at
its default value of $0.02$ seconds.

PhysX uses the Coulomb model for friction and restitution is a measure of the
loss of kinetic energy between colliding particles. Details and further
references about the rigid body system and friction model used in PhysX can be
found in the work by Tonge [16].

The model incorporates the following parameters which are freely specifiable
and chosen to agree with experimental results: friction, restitution, and
sphere diameter. In addition, our model allows us to vary the angle of
inclination of the plane, the gap width and the distance from the obstruction
to the gap. We incorporate imperfections in the sphere diameters by allowing
them to vary by up to 5%. The model uses $1551$ particles whilst $\approx
1000$ particles were used in the experiment. Each simulation is performed $20$
times.

Unity [13], a freely available game engine, is used as our programming
interface to PhysX. We also use Blender [14], a freely available 3D modeling
package for modeling the obstructions.

##  3 Flow rate without obstructions

The basic set up of the experiment with obstructions is shown in Fig. 1(a),
where particles are released from above the obstruction. The particles
themselves are spheres with diameters: $6.15$ mm, $7.75$ mm and $11.9$ mm,
referred to as ‘small’, ‘medium’, and ‘large’ respectively, and are shown in
Fig. 1 (b). Note that without obstructions, the particles are released at the
gap as shown in Fig. 1(c).

Figure 1: (a) Basic set up of the experiment with an obstruction. (b)
Particles are spheres of three different diameters. (c) Experiment of flow
rate without obstructions, the particles are released at the gap.

We use an inclined plane that makes an angle of $8^{o}$ with the horizontal.
The plane has retaining walls, walls that form a bottleneck and release
barriers as shown in Fig. 1(a). We pack a known number of spheres of a
specified size, shown without texture in Fig. 1(b), onto the plane and release
them from rest so that they collide with the bottleneck walls, and flow
through the gap. The walls of the bottleneck make an angle of $30^{o}$ to the
retaining walls.

Under gravity, the spheres move down the inclined plane, come into contact
with the bottleneck walls and flow through the gap. The number of particles is
known and the time taken for all particles to flow through the gap is
measured, from which the flow rate, $J$, is determined. The flow rate is thus
defined as the number of particles passing through the gap per second.

The experiment is repeated for various gap widths, with the corresponding flow
rates measured. The entire experiment is performed individually for small,
medium, and large spheres. All dimensions (particles, plane, bottleneck
angles, etc) are known, as are the coefficients of friction and restitution of
the particles. As per the experiment, in our model we release the particles at
rest from the gap as shown in Fig. 1 (c).

It was found that there is a relatively sharp rise in flow rate as the gap is
increased. It is precisely this range of distances where the flow rate
transitions from clogging to intermittent to continuous flow rates. The model
also exhibits clogging as shown in Fig. 2 (a), and intermittent flow.
Therefore, in order to get reliable and reproducible flow rates, we decided to
perturb the particles just sufficiently to prevent extended periods of
clogging which thereby affect flow rates. As described by Garcimartin et al.
[11], the problem of clogging can effectively be eliminated by applying
external vibrations.

Figure 2: The model exhibits clogging (a) at the gap and between the
obstruction and the bottleneck walls (b), consistent with the experiment.
Figure 3: (a) directions of shaking of the bottleneck walls in the experiment
of flow rate without obstructions. (b) directions of shaking in the experiment
of flow rate with obstructions. Shaking is necessary to avoid clogging between
the obstruction and the bottleneck walls. (c) obstructions used in
experiments. The V obstruction is the A obstruction rotated through $180^{o}$.

In a similar way, we use the technique of ‘shaking’. In the case of flow rate
without obstructions, the bottleneck walls were shaken as shown in Fig. 3(a)
from a maximum of $0.5$ mm diminishing to zero over one second, every five
seconds. The parameters that define a shake are: shake amplitude, direction of
shake, shake duration and time between shakes. These parameters can be set
arbitrarily in the model and, whilst other times and distances were tested,
the above were found to be just sufficient to yield reproducible results.

We varied the values of coefficient of friction and coefficient of restitution
in the model until we got good agreement with experimental results. It was
found that varying the coefficient of restitution had little to no effect on
flow rate.

We found a coefficient of friction of $0.24$ yielded model results that gave
good agreement with experimental results. The experimentally measured value of
the spheres was found to be $0.28\pm 0.4$. The results for each diameter
sphere are shown in Fig. 4. The model results agree with the experimental
results quite well.

Figure 4: Flow rate without obstruction versus gap in dimensionless form.
Comparison of experiment with simulation using spheres of three different
diameters. $J$ is the flow rate, $g$ is the acceleration due to gravity, $d$
is the particle diameter and $W$ is the gap width.

##  4 Flow rate with obstructions

The experiment is repeated with medium spheres (of diameter $7.75$ mm).
However, this time, obstructions (shown in Fig. 3 (c)) are placed at varying
distances away from the gap, which is now fixed at a width of $3.3$ cm.

It was found that there was a relatively sharp rise in flow rate as an
obstruction is moved away from the gap. It is precisely this range of
distances where the flow rate transitions from clogging to intermittent to
continuous flow rates. The model also exhibited clogging as shown in Fig.
2(b), and intermittent flow. Therefore, in order to get reliable and
reproducible flow rates, we decided to perturb the particles just sufficiently
to prevent extended periods of clogging which thereby affect flow rates.

With the gap set at $3.3$ cm, no clogging was observed for spheres of $7.75$
mm at the gap. However, with the presence of an obstruction, clogging was
observed in areas between the obstruction and the bottleneck walls as shown in
Fig. 2(b) and Fig. 3(b). To get consistently measurable flow rates in the
model, once again, we decided to perturb the system using the technique of
‘shaking’, described above. This time, we decided to shake the obstruction
also shown in Fig. 3(b). The obstruction was shaken by a maximum of $1$ mm
diminishing to zero over one second, every five seconds. The freely
specifiable parameters that define a ‘shake’ are the same as those for shaking
the bottleneck walls. Figure 5 shows the maximum flow rates for the three
obstructions, where the ‘waiting room’ effect, discussed in Ref. [8], occurs
for the cylinder and the V obstruction but not the A obstruction.

Figure 5: Peak flow rates with obstructions, where the ‘waiting room’ effect
is present for the cylinder and the V obstruction but not for the A
obstruction. Figure 6: Results of flow rate with obstructions versus distance
in dimensionless form. $J$ is the flow rate, $g$ is the acceleration due to
gravity, $d$ is the particle diameter, and $D$ is the distance from the
obstruction to the gap. The experimental results in (a) show that each
obstruction has a range of distances where flow rate is greater than with no
obstruction. The numerical results (b-d) show that the model also predicts
that each obstruction has a range of distances from the gap where the flow
rate is greater than flow rate with no obstruction.

Figure 6 shows flow rates as measured experimentally and with the model. Both
sets of results indicate increases in flow rates for certain ranges of
distance of obstruction to the gap. For the model, each graph shows error bars
and data that were more than two standard deviations away from the mean were
not included in the statistical analysis. This was necessary as, even with the
shaking technique described above, the particles did occasionally clog for
extended periods of time, thereby affecting the flow rate measurements.

Clearly, the model shows all obstructions improve flow rate for a range of
distances, and reassuringly, the further the obstruction is moved away from
the gap, the closer we get to flow rates without obstruction.

In the case of the cylinder, due to its shape and size relative to the
bottleneck, we could only start at distances $\sim 7$ cm from the gap. We see
that there is a steady rise in flow rate as the cylinder is moved further away
from the gap. A maximum flow rate is reached and then a slightly sharper
decrease is reached until we get flow rates the same as if no obstruction is
present.

In the case of the A obstruction, due to its shape and size relative to the
bottleneck, we were able to place the obstruction very close to the gap, which
is why distances start at $\sim 0.4$ cm from the gap. We see a sharp rise in
flow rate as the obstruction is moved away from the gap, and then we see a
leveling out of the flow rate over a range of distances from $0.8$ cm to $2.8$
cm. There is a small peak in this range, but we consider this a statistical
fluctuation.

Moving further away, from $2.8$ cm to $3.2$ cm, we see a sharp decrease in
flow rate to values that are actually below flow rates with no obstacle. This
is an intriguing phenomenon, showing what we might expect: obstructions
decrease flow rates.

The flow rate remains constant to a distance of $3.6$ cm and then starts
rising. At a distance of $4$ cm, the flow rate is approximately equal to that
with no obstructions.

In the case of the V obstruction, due to its shape and size relative to the
bottleneck, we could only start at distances $\sim 6$ cm from the gap. We see
that there is a steady rise in flow rate that is slower than that of the A
obstruction, as the V obstruction is moved further away from the gap. A
maximum flow rate is reached, higher than that with no obstruction, and then a
decrease occurs until we get flow rates the same as if no obstruction were
present.

##  5 Limitations of PhysX engine

We observed that the model exhibits continuous flow, intermittent flow, and
clogging, all phenomena that have been identified and observed experimentally.
The model was quite sensitive to the absolute sizes of the geometric
structures used. This is true for all physics engines, as they have to be
‘tuned’ to the range we would like them to work most accurately in. For
example, an engine might be tuned to work best with objects whose sizes are of
the order of meters but it will not work for objects at the scale of
millimeters. To overcome this limitation, we ran the simulations in the
dimensions where the performance of the engine was optimal, and then used
dimensional analysis to extrapolate the result to the real experimental
conditions.

We found that the particles in the model needed to be perturbed more often
than in actual experiments, in order to get consistent flow rates. This can be
attributed to a real time optimization employed by PhysX, known as ‘sleeping’.
When a particle’s angular and/or linear velocity falls below certain threshold
values for more than a few frames, these velocities are set to zero and the
particle goes into a ‘sleep’ state, in which no collision detection occurs and
hence the particle’s velocity remains at zero. The particle ‘wakes up’ when it
is subjected to net forces. The reason for this optimization is that it
relieves the processor of having to perform needless computations especially
with regard to collision detection, thereby allowing much larger numbers of
particles to be present and only performing the necessary computations as
required. For this study, we left the threshold ‘sleep velocities’ at their
default values.

In spite of these limitations, we can confirm that, as reported in various
experiments, it is possible to increase flow rate of discrete particles by
placing an obstruction at a suitable location near a gap.

##  6 Conclusions

We have developed a model that exhibited reasonable quantitative accuracy in
the case of no obstructions, and good qualitative agreement when obstructions
were introduced. As in experiments, we showed that the flow exhibits clogging,
intermittent flow rates, and continuous flow rates. The model also confirms
experimental results that placement of an obstruction near a gap can actually
increase flow rate through the gap. Better experimental agreement with
obstructions will be possible if the roundness of the obstructions is more
accurately modeled.

The model is quite flexible in that many parameters can be changed via user
input, and without modifying the actual code. In this way, the model can be
used as a convenient tool to suggest further experiments and allows us to
investigate different obstruction shapes, both of which may help us gain a
greater understanding of granular flow.

We have also described a suite of freely available realtime 3D software, which
is mature and has many online resources in the form of documentation and
tutorials and is being actively developed, which can be used together to
create simulations of interest in the area of granular flow.

###### Acknowledgements.

We thank Ivan Gojkovic and Jaeris Wu for performing the experimental part of
this paper.

## References

  * [1] C Magalhaes, A Atman, G Combe, J Moreira,  Jamming transition in a two-dimensional open granular pile with rolling resistance, Pap. Phys. 6, 060007 (2014).
  * [2] W Beverloo, H Leniger, J van der Velde,  The flow of granular solids through orifices, Chem. Eng. Sci. 15, 260 (1961).
  * [3] T Wilson, C Pfeifer, N Mesyngier, D. Durian,  Granular discharge rate for submerged hoppers, Pap. Phys. 6, 060009 (2014).
  * [4] I Zuriguel,  Clogging of granular materials in bottlenecks, Pap. Phys. 6, 060014 (2014).
  * [5] P Mort,  Characterizing flowability of granular materials by onset of jamming in orifice flows, Pap. Phys. 7, 070004 (2015).
  * [6] O Rodolfo, J Sales, M Gargiulo, A Vidales,  Density distribution of particles upon jamming after an avalanche in a 2D silo, Pap. Phys. 7, 070007 (2015).
  * [7] G Lumay, J Schockmel, D Hernandez-Enriquez, S Dorbolo, N Vandewalle, F Pacheco-Vazquez,  Flow of magnetic repelling grains in a two-dimensional silo, Pap. Phys. 7, 070013 (2015).
  * [8] F Alonso-Marroquin, S Azeezullah, S Gallindo-Torres, L Olsen-Kettle,  Bottlenecks in granular flow: When does an obstacle increase flow rate in an hourglass?, Phys. Rev. E 85, 020301 (2012).
  * [9] S Yang, S Hsiau,  The simulation and experimental study of granular materials discharged from a silo with the placement of inserts, Powder Tech. 120, 244 (2001).
  * [10] I Zuriguel, A Janda, A Garcimartín, C Lozano, R Arévalo, D Maza, Silo clogging reduction by the presence of an obstacle,  Phys. Rev. Lett. 107, 278001 (2011).
  * [11] A Garcimartín, C Lozano, G Lumay, I Zuriguel,  Avoiding clogs: The shape of arches and their stability against vibrations,  AIP Conf. Proc. 1542, 686 (2013).
  * [12] K Kumar,  Learning physics modeling with PhysX,  Packt Publishing Ltd, Birmingham (2013).

  * [13] W Goldstone,  Unity game development essentials,  Packt Publishing Ltd, Birmingham (2009).
  * [14] J Blain,  The complete guide to Blender graphics, (Second Edition),  A K Peters/CRC Press, Boca Raton (2014).
  * [15] C Carlevaro, L Pugnaloni,  Steady state of trapped granular polygons,  J. Stat. Mech. P01007 (2011).
  * [16] R Tonge,  Iterative rigid body solvers,  `http://www.gdcvault.com/play/1018160/`
`Physics-for-Game`.


# A Dynamic Data Structure for Representing Temporal Transitive Closures on
Disk

Luiz F. Afra Brito  Federal University of Uberlândia, Brazil  Marcelo
Albertini  Federal University of Uberlândia, Brazil  Bruno A. N. Travençolo
Federal University of Uberlândia, Brazil

###### Abstract

Temporal graphs represent interactions between entities over time. These
interactions may be direct, a contact between two vertices at some time
instant, or indirect, through sequences of contacts called journeys. Deciding
whether an entity can reach another through a journey is useful for various
applications in complex networks. In this paper, we present a disk-based data
structure that maintains temporal reachability information under the addition
of new contacts in a non-chronological order. It represents the _timed
transitive closure_ (TTC) by a set of _expanded_ R-tuples of the form
$(u,v,t^{-},t^{+})$, which encodes the existence of journeys from vertex $u$
to vertex $v$ with departure at time $t^{-}$ and arrival at time $t^{+}$. Let
$n$ be the number of vertices and $\tau$ be the number of timestamps in the
lifetime of the temporal graph. Our data structure explicitly maintains this
information in linear arrays using $O(n^{2}\tau)$ space so that sequential
accesses on disk are prioritized. Furthermore, it adds a new unsorted contact
$(u,v,t)$ accessing $O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$ sequential
pages in the worst-case, where $B$ is the of pages on disk; it answers whether
there is of a journey from a vertex $u$ to a vertex $v$ within a time interval
$[t_{1},t_{2}]$ accessing a single page; it answers whether all vertices can
reach each other in $[t_{1},t_{2}]$; and it reconstructs a valid journey that
validates the reachability from a vertex $u$ to a vertex $v$ within
$[t_{1},t_{1}]$ accessing $O\left(\nicefrac{{n\tau}}{{B}}\right)$ pages. Our
experiments show that our novel data structure are better that the best known
approach for the majority of cases using synthetic and real world datasets.

##  1 Introduction

Temporal graphs represent interactions between entities over time. These
interactions often appear as contacts at specific timestamps. Entities can
also interact indirectly with each other by chaining several contacts. For
example, in a communication network, devices that are physically connected can
send new messages or propagate received ones; thus, by first sending a new
message and repeatedly propagating messages over time, remote entities can
communicate indirectly. Time-respecting paths in temporal graphs are known as
temporal paths, or simply journeys, and when a journey exists from one vertex
to another one, we say that the first can reach the second.

In a computational environment, it is often useful to check whether entities
can reach each other [19, 6, 21, 23, 22, 3, 15, 5]. Beyond the sole
reachability, some applications also require the ability to reconstruct a
journey if one exists [23, 11, 25, 13, 14, 16]. In standard graphs, the
problem of updating reachability information is known as _dynamic
connectivity_ and it has been extensively studied [17, 12, 9, 26, 18, 20]. In
temporal graphs, fewer studies addressed this problem. For instance, in [2,
21], the authors assume that input is chronologically ordered and they give
worst-case complexities. In [24], the authors assume non-chronological input,
whereas their strategy is optimized for the average case.

Particularly to our interest, in [4], the authors considered the _only-
incremental_ problem, which supports only the addition of unsorted contacts.
Their data structure supports the following four operations, where, by
convention, $\mathcal{G}$ is a temporal graph, $u$ and $v$ are vertices of
$\mathcal{G}$, and $t,t_{1}$, and $t_{2}$ are timestamps:

  1. 1.

add_contact(u,v,t), which updates information based on a contact from $u$ to
$v$ at time $t$;

  2. 2.

can_reach($u$,$v$,$t_{1}$,$t_{2}$), which returns true if $u$ can reach $v$
within the interval $[t_{1},t_{2}]$;

  3. 3.

is_connected($t_{1}$,$t_{2}$), which returns true if $\mathcal{G}$ restricted
to the interval $[t_{1},t_{2}]$ is temporally connected, i.e., all vertices
can reach each other within the interval $[t_{1},t_{2}]$; and

  4. 4.

reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$), which returns a journey (if one
exists) from $u$ to $v$ occurring within the interval $[t_{1},t_{2}]$.

Their update algorithm maintains a timed transitive closure (TTC), a concept
that generalizes the transitive closure for temporal graphs based on
reachability tuples (R-tuples), in the form $(u,v,t^{-},t^{+})$, representing
journeys from vertex $u$ to $v$ departing at $t^{-}$ and arriving at $t^{+}$.
Their data structure uses $O\left(n^{2}\tau\right)$ space while supporting
add_contact, can_reach, is_connected, and reconstruct_journey, respectively,
in $O\left(n^{2}\log\tau\right)$, $O\left(\log\tau\right)$,
$O\left(n^{2}\log\tau\right)$, and $O\left(k\log\tau\right)$ worst-case time,
where $n$ is the number of vertices of the temporal graph, $\tau$ is the
number of time instances, and $k$ is the length of the resulting journey.

However, they keep their data structure in primary memory and the cost of
storing and maintaining large TTCs is prohibitive. We conducted a simple
experiment to show how much space is necessary for temporal reachability.
First, we generated random temporal graphs using the Edge-Markovian Evolving
Graph (EMEG) model [8]. In this model, if an edge is active at time $t-1$,
then it has probability $p$ of disappearing at time $t$, otherwise, it has
probability $q$ of appearing at time $t$. We represented temporal graphs in
memory using adjacency matrices storing, in each cell, timestamps at which
edges are active. Then, we built the corresponding TTCs using the approach
described in [4]. In this experiment, we varied the number of vertices $n$ and
the number of time instances $\tau$ while fixing $p=0.1$ and $q=0.3$.

In Table 1, we see, for example, that a temporal graph with $512$ vertices and
$\tau=64$ produced by the EMEG model has $2.8$ million contacts, and we needed
around $33$ MBs of space to store it in memory. Besides, we needed around
$156$ MBs of space to store the corresponding TTC, which, in this case, it is
almost five times the space needed to store the temporal graph.

$n$ | $\tau$ | $|C|$ | $data(\mathcal{G})$ | $data({\text{\footnotesize TTC}})$
---|---|---|---|---
32 | 8 | 1268 | 0.02 | 0.01
32 | 16 | 2670 | 0.03 | 0.12
32 | 32 | 5249 | 0.06 | 0.24
64 | 8 | 5539 | 0.08 | 0.27
64 | 16 | 10908 | 0.14 | 0.54
64 | 32 | 21421 | 0.26 | 1.08
64 | 64 | 42671 | 0.50 | 2.17
128 | 8 | 21203 | 0.31 | 1.13
128 | 16 | 43011 | 0.55 | 2.30
128 | 32 | 86746 | 1.06 | 4.63
128 | 64 | 173479 | 2.05 | 9.31
256 | 8 | 86574 | 1.24 | 4.67
256 | 16 | 174970 | 2.25 | 9.51
256 | 32 | 346994 | 4.22 | 19.17
256 | 64 | 696436 | 8.22 | 38.54
512 | 8 | 349114 | 5.00 | 19.01
512 | 16 | 702294 | 9.04 | 38.66
512 | 32 | 1396033 | 16.98 | 78.00
512 | 64 | 2800520 | 33.05 | 156.64
Table 1: Space for storing temporal graphs with $n$ vertices, $\tau$ time
instances and $|C|$ contacts, and their corresponding TTCs. Columns
$data(\mathcal{G})$ and $data({\text{\footnotesize TTC}})$ represent,
respectively, the space in megabytes of the generated temporal graphs and
their TTCs.

Next, we built a linear regression model with the data presented in Table 1 in
order to extrapolate the input parameters. Consider, for example, the scenario
in which one million people use a bluetooth device that registers when and who
gets close to each other and sends this information to a centralized server.
Consider also that each individual makes in average $30$ contacts per day. In
this setting, by using our model, we could check that a centralized server
would require at least $100$ GBs of space in less than a year to store just
the plain contacts as a temporal graph. If one needs to support reachability
queries by using a TTC, it would be necessary roughly $600$ GBs of space.

Motivated by such scenarios, we investigate the problem of maintaining TTCs on
disk. A simple, but not efficient, approach would be to naively implement a
data structure for disk based on the approach described in [4]. Briefly, their
strategy maintains self-balanced binary search trees (BSTs) containing time
intervals for each pair of vertices in order to retrieve reachability
information. However, this approach does not consider data locality, thus each
update operation would randomly access an excessive amount of pages on disk to
retrieve information from each BST. For instance, if we use B+-trees [1] as a
replacement for BSTs, their algorithm for answering add_contact would access
$O\left(n^{2}\right)$ B+-trees and, in each B+-tree, it would access
$O\left(\log_{B}{\tau}\right)$ pages, where $B$ is the page size, resulting in
$O\left(n^{2}\log_{B}{\tau}\right)$ random accesses on disk. Therefore, we
need a novel approach that better organizes data on disk.

We propose in this paper an incremental disk-based data structure that reduces
the number of disk accesses for both update and query operations while
prioritizing sequential accesses. The core idea of our novel approach is to
maintain explicitly an expanded set of non-redundant R-tuples containing
$n^{2}\tau$ elements. Conceptually, we maintain it using two $3$-dimensional
arrays, $M_{out}$ and $M_{in}$, of size $n\times\tau\times n$ such that
$M_{out}[u,t^{-},v]=t^{+}$ and $M_{in}[v,t^{+},u]=t^{-}$. The former supports
querying the earliest arrival time $t^{+}$ a journey departing from vertex $u$
at time $t^{-}$ can arrive at vertex $v$, and the latter supports querying the
latest departure time $t^{-}$ a journey arriving to vertex $v$ at time $t^{+}$
can depart from vertex $u$.

Our algorithm to compute add_contact eagerly updates both arrays accessing
$O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$ disk pages in the worst case.
Despite having a linear factor on $\tau$ instead of logarithmic, the expected
cost of our update routine reduces considerably as we insert new contacts.
This is because journey schedules become stricter and the probability of
replacing them with faster ones reduces. Since we explicitly maintain
reachability information, our algorithms to answer can_reach, is_connected,
and reconstruct_journey access, respectively, one,
$\Theta\left(\nicefrac{{n^{2}}}{{B}}\right)$ and
$\Theta\left(\nicefrac{{n}}{{B}}\right)$ pages.

We compare our novel data structure with a naïve adaptation of the approach
introduced in [4] using B+-trees as replacement for BSTs. Our experiments show
that our novel data structure performs better on the synthetic datasets and on
the majority of real-world datasets we used. Even though the worst-case
complexity of our algorithm for the add_contact(u,v,t)uvt operation is linear
in $\tau$ instead of logarithmic, it runs much faster on average. We attribute
this behavior to the fact that as new contacts are inserted, our data
structure updates on average only a few cells of both arrays $M_{out}$ and
$M_{in}$.

We organized this paper as follows. In Section 2, we present the definitions
used throughout this paper. In Section 3, we define our expanded set of
R-tuples, introduce our new data structure to represent TTCs on disk, and
provide low-level primitives for manipulating them. In Section 4, we describe
our algorithms for each operation using our data structure along with their
complexities in terms of number of disk accesses. In Section 5, we investigate
the execution of our algorithms by comparing them with our implementation
using B+-trees. Finally, Section 6 concludes with some remarks and open
questions.

##  2 Definitions

Following the definition in [7], a temporal graph is a tuple
$\mathcal{G}=(V,E,\mathcal{T},\rho,\zeta)$. Sets $V$ and $E\subseteq{}V\times
V$ represent the vertices and the edges of the underlying standard graph.
Interval $\mathcal{T}=[1,\tau]\subset\mathbb{N}$ describes the lifetime of the
temporal graph. We consider in this paper that $E$ is a set of directed edges.
Functions $\rho:E\times\mathcal{T}\to\\{0,1\\}$ and
$\zeta:E\times\mathcal{T}\mapsto\mathbb{N}$ are, respectively, the _presence
function_ and the latency function. The presence function expresses whether an
edge is present at a time instant. We also call $(u,v,t)$ a _contact_ in
${\mathcal{G}}$ if $\rho((u,v),t)=1$. The latency function expresses the
duration of an interaction for an edge at a time. Here, we use a constant
latency function, namely $\zeta=\delta$, where $\delta$ is any fixed positive
integer.

We define reachability in temporal graphs in a time-respecting way, by
requiring that a path travels along non-decreasing ($\delta=0$) or increasing
($\delta\geq 1$) times. These paths are called temporal paths or journeys
interchangeably.

######  Definition 1 (Journey).

A journey from $u$ to $v$ in $\mathcal{G}$ is a sequence of contacts
$\mathcal{J}=\langle c_{1},c_{2},\ldots,c_{k}\rangle$, whose sequence of
underlying edges form a valid time-respecting path from $u$ to $v$. For each
contact $c_{i}=(u_{i},v_{i},t_{i})$, it holds that
$\rho((u_{i},v_{i}),t_{i})=1$, $v_{i}=u_{i+1}$, and $t_{i+1}\geq t_{i}+\delta$
for $i\in[1,k-1]$. We say that $departure(\mathcal{J})=t_{1}$,
$arrival(\mathcal{J})=t_{k}+\delta$ and
$duration(\mathcal{J})=arrival(\mathcal{J})-departure(\mathcal{J})$. A journey
is _trivial_ if it comprises a single contact.

######  Definition 2 (Reachability).

A vertex $u$ can _reach_ a vertex $v$ within time interval $[t_{1},t_{2}]$ iff
there is a journey $\mathcal{J}$ from $u$ to $v$ in $\mathcal{G}$ that departs
at $departure(\mathcal{J})\geq t_{1}$ and arrives at $arrival(\mathcal{J})\leq
t_{2}$.

Just as the number of paths in a standard graph, the number of journeys in a
temporal graph could be too large to be stored explicitly (typically,
factorial in $n$). To avoid this problem, R-tuples capture the fact that a
vertex can reach another one within a certain time interval without storing
the corresponding journeys [4].

######  Definition 3 (R-tuple).

A R-tuple is a tuple $r=(u,v,t^{-},t^{+})$, where $u$ and $v$ are vertices in
$\mathcal{G}$, and $t^{-}$ and $t^{+}$ are timestamps in $\mathcal{T}$. It
encodes the fact that vertex $u$ can reach vertex $v$ through a journey
$\mathcal{J}$ such that $departure(\mathcal{J})=t^{-}$ and
$arrival(\mathcal{J})=t^{+}$. If several such journeys exist, then they are
all represented by the same R-tuple.

Lastly, given a temporal graph $\mathcal{G}$, the timed transitive closure
(TTC) of $\mathcal{G}$ is a directed multigraph on the same set of vertices,
whose edges correspond to the minimal, i.e., non-redundant, set of R-tuples of
$\mathcal{G}$. The purpose of TTCs is to encode reachability information among
vertices, parametrized by time intervals, so that one can subsetuently decide
if a new contact can be composed with existing journeys. In paper [4], the
authors showed that there are $O\left(n^{2}\tau\right)$ non-redundant R-tuples
in a temporal graph $\mathcal{G}$ and it comprises those whose intervals do
not include each other for the same pair of vertices, i.e., only information
regarding the fastest journeys.

##  3 Disk-Based Timed Transitive Closure

In this section, we describe our novel approach to maintain TTCs in secondary
memory. First, in Section 3.1, we define the concept of an _expanded_ set of
representative R-tuples and show that it has size
$\Theta\left(n^{2}\tau\right)$. Then, in Section 3.2, we introduce our new
data structure that uses this expanded set in order to improve the maintenance
of data in non-uniform access storages and provide direct access to
reachability information.

###  3.1 Expanded Reachability Tuples (Expanded R-tuples)

The data structure introduced in [4] spreads the minimal set of R-tuples into
multiple BSTs, each one concerning a unique pair of vertices. The authors
store these BSTs in separated regions of memory and, therefore, the
organization of data is not optimal when working with storages that have non-
uniform access time.

In order to mitigate this problem, we define an expanded set of R-tuples
$(u,v,t^{-},t^{+})$ that is easier to maintain sequentially, since we can use
continuous arrays indexed by $t^{-}$ or $t^{+}$. First, we define the _left_
and _right_ expansion of a single R-tuple.

######  Definition 4 (Left and right expansion).

The left expansion of a R-tuple $r=(u,v,t^{-},t^{+})$ is the set containing
all R-tuples $(u,v,t,t^{+})$ for $1\leq t\leq t^{-}$. Similarly, the right
expansion of $r$ is the set containing all R-tuples $(u,v,t^{-},t)$ for
$t^{+}\leq t\leq\tau+\delta$.

The R-tuples produced by the left expansion of a R-tuple $r$ are valid because
a source vertex departing earlier can simply wait until the departure time of
$r$, and take the original journey described by $r$. Similarly, the R-tuples
produced by the right expansion of $r$ are valid because, after taking the
original journey described by $r$, a destination vertex can simply wait until
the arrival time of the new R-tuple.

Applying both expansions to each R-tuple in a set $\mathcal{R}$ and taking the
union of the sets produced by the same expansion creates two separated
expanded sets, the _left-expanded_ set $\mathcal{R}_{\text{\scriptsize
left}}$, and the _right-expanded_ set $\mathcal{R}_{\text{\scriptsize
right}}$. For each expanded set, we define an inclusion operator.

######  Definition 5 (Left and right inclusion).

Given any two R-tuples $r_{1}=(u_{1},v_{1},t^{-}_{1},t^{+}_{1})$ and
$r_{2}=(u_{2},v_{2},t^{-}_{2},t^{+}_{2})$ in $\mathcal{R}_{\text{\scriptsize
left}}$, $r_{1}\subseteq_{\text{\scriptsize left}}r_{2}$ if and only if
$u_{1}=u_{2}$, $v_{1}=v_{2}$, $t^{-}_{1}=t^{-}_{2}$, and $t^{+}_{1}\leq
t^{+}_{2}$. Similarly, if $r_{1}$ and $r_{2}$ are in
$\mathcal{R}_{\text{\scriptsize right}}$, $r_{1}\subseteq_{\text{\scriptsize
right}}r_{2}$ if and only if $u_{1}=u_{2}$, $v_{1}=v_{2}$, $t^{-}_{1}\geq
t^{-}_{2}$, and $t^{+}_{1}=t^{+}_{2}$.

However, R-tuples produced by expansion can share redundant information. For
example, consider the R-tuples $r_{1}=(a,b,2,7)$ and $r_{2}=(a,b,2,9)$. Both
R-tuples represent journeys that departs from vertex $a$ at time $2$ and
arrives at vertex $b$, one at time $7$ and the other at time $9$. In this
case, $r_{2}$ can be safely discarded since we can take a journey represented
by $r_{1}$ ending at time $7$ and wait at vertex $v$ until time $9$.
Redundancy of R-tuples in $\mathcal{R}_{\text{\scriptsize left}}$ and
$\mathcal{R}_{\text{\scriptsize right}}$ is treated differently using their
corresponding inclusion operators.

######  Definition 6 (Left and right redundancy).

Let $r\in\mathcal{R}_{\text{\scriptsize left}}$, $r$ is called left-redundant
in $\mathcal{R}_{\text{\scriptsize left}}$ if there is
$r^{\prime}\in\mathcal{R}_{\text{\scriptsize left}}$ such that
$r^{\prime}\subseteq_{\text{\scriptsize left}}r$. Similarly, if
$r\in\mathcal{R}_{\text{\scriptsize right}}$, $r$ is called right-redundant in
$\mathcal{R}_{\text{\scriptsize right}}$ if there is
$r^{\prime}\in\mathcal{R}_{\text{\scriptsize right}}$ such that
$r^{\prime}\subseteq_{\text{\scriptsize right}}r$. A set
$\mathcal{R}^{*}_{\text{\scriptsize left}}$ with no left-redundant R-tuple is
called left non-redundant and a set $\mathcal{R}^{*}_{\text{\scriptsize
right}}$ with no right-redundant R-tuple is called right no-redundant.

######  Lemma 1.

The maximum size of a left non-redundant or right non-redundant set of
R-tuples for $\mathcal{G}$ is $\Theta\left(n^{2}\tau\right)$.

###### Proof.

It suffices to prove that the maximum number of pairwise incomparable R-tuples
in $\mathcal{R}_{\text{\scriptsize left}}$ and $\mathcal{R}_{\text{\scriptsize
right}}$ is $O\left(n^{2}\tau\right)$, since some graphs induce
$\Theta\left(n^{2}\tau\right)$ incomparable R-tuples from unexpanded sets, see
[4]. We prove that the sizes of $\mathcal{R}_{\text{\scriptsize left}}$ and
$\mathcal{R}_{\text{\scriptsize right}}$ are $O\left(n^{2}\tau\right)$ as
follows. There are $\Theta\left(n^{2}\right)$ ordered pairs of vertices. Thus,
it is enough to show that for each pair $(u,v)$, the number of incomparable
R-tuples in $\mathcal{R}^{u,v}_{\text{\scriptsize left}}$ and
$\mathcal{R}^{u,v}_{\text{\scriptsize right}}$, whose source vertex is $u$ and
destination vertex is $v$, is $\Theta\left(\tau\right)$. Let
$\mathcal{R}^{u,v}_{\text{\scriptsize left}}$ be a left non-redundant set of
such R-tuples, as every incomparable R-tuple has different arrival timestamps,
$|\mathcal{R}_{\text{\scriptsize left}}|\leq\tau$. Similarly, let
$\mathcal{R}^{u,v}_{\text{\scriptsize right}}$ be a right non-redundant set of
such R-tuples, as every incomparable R-tuple has different departure
timestamps, $|\mathcal{R}_{\text{\scriptsize right}}|\leq\tau$.
∎

###  3.2 Encoding the TTC on Disk

|  |  |
---|---|---|---
(a) $\mathcal{G}$ |  | (b) $M_{out}$ | (c) $M_{in}$
Figure 1:  Temporal graph and its associated reachability data structure. In
(a), we show a temporal graph with three vertices. Numbers on edges represent
the time in which edges are active. Edges with the same color form a journey
from vertex $u$ to vertex $v$. In (b), we show the corresponding arrays
$M_{out}$ and $M_{in}$ considering $\delta=1$. Both arrays are depicted as
2-dimensional arrays by grouping their first two dimensions. For instance,
$M_{out}[u,2,w]=M_{out}[(u,2),w]=3$. Cells have the same color as the
contacts, i.e., the edge at a timestamp, that originated the update. $M_{out}$
stores the minimum possible arrival timestamps to destinations and $M_{in}$
sotores the maximum possible departure timestamps from origins.

We encode the TTC using two $3$-dimensional arrays, $M_{out}[u,t^{-},v]=t^{+}$
and $M_{in}[v,t^{+},u]=t^{-}$, both with dimensions $n\times\tau\times n$,
representing expanded sets of R-tuples. Each cell in $M_{out}$ represents a
R-tuple in $\mathcal{R}^{*}_{\text{\scriptsize left}}$ by storing the earliest
possible arrival time $t^{+}$ at which a vertex $u$ departing at time $t^{-}$
can reach a vertex $v$ through a journey. If there is a cell
$M_{out}[u,t^{-},v]=t^{+}$, then all cells $M_{out}[u,t,v]$, for
$t\in[1,t^{-}-1]$ must have an arrival time $t_{\text{\scriptsize left}}\leq
t^{+}$, since a journey from $u$ departing at a time $t<t^{-}$ can simply wait
at vertex $u$ until time $t^{-}$ and then use the remaining path already
described by $M_{out}[u,t^{-},v]=t^{+}$. Similarly, each cell in $M_{in}$
represents a R-tuple in $\mathcal{R}^{*}_{\text{\scriptsize right}}$ by
storing the latest possible departure time $t^{-}$ at which a vertex $v$ can
arrive at time $t^{+}$ to a vertex $u$ through a journey. If there is a cell
$M_{in}[v,t^{+},u]=t^{-}$, then all cells $M_{in}[v,t,u]$, for
$t\in[t^{+}+1,\tau+\delta]$, must have a departure time $t_{\text{\scriptsize
right}}\geq t^{+}$, since a journey to $v$ arriving at a time $t>t^{+}$ can
use the path already described by $M_{in}[v,t^{+},u]=t^{-}$ and then simply
wait at vertex $v$ until time $t$. During the creation of a TTC, $M_{out}$
cells are initialized with $\infty$ and $M_{in}$ cells with $-\infty$. Figure
1 illustrates both $M_{out}$ and $M_{in}$.

Internally, we represent $M_{out}$ and $M_{in}$ as one-dimensional arrays
using, respectively, the mapping functions $F_{out}\colon(u,t^{-},v)\mapsto
n(u\tau+\tau-(t^{-}+1))+v+1$ and $F_{in}\colon(v,t^{+},u)\mapsto
n(v\tau+t^{+}-\delta)+u+1$. Observing Figure 1, $F_{out}$ arranges the cells
of $M_{out}$ by row (left to right) and, for each source vertex, later
departures come first. $F_{in}$ also arranges $M_{in}$ by row but, in
contrast, for each destination vertex, earlier arrivals come first. By
subtracting $\delta$ from $t^{+}$ in $F_{in}$, we ensure all $t^{+}$ values
fit in $M_{in}$. Thus, reading sequentially the range
$[F_{out}(u,t^{-},1),F_{out}(u,t^{-},n)]$ from $M_{out}$ gives direct access
to the earliest arrival times to reach all vertices when departing from $u$ at
time $t^{-}$. Similarly, reading sequentially the range
$[F_{in}(v,t^{+},1),F_{in}(v,t^{+},n)]$ from $M_{in}$ gives direct access to
the latest departure times to leave all vertices when arriving at $v$ at time
$t^{+}$.

Finally, assuming a general function $F$ that maps to $F_{out}$, whether
accessing $M_{out}$, or $F_{in}$, whether accessing $M_{in}$, we provide the
following low-level operations for manipulating our data structures on disk:

  1. 1.

read_cell($M,w_{1},t,w_{2}$), which returns the value of $M$ ($M_{out}$ or
$M_{in}$) at position $F(w_{1},t,w_{2})$;

  2. 2.

write_cell($M,w_{1},t_{1},w_{2},t_{2}$), which replaces the value of $M$ at
position $F(w_{1},t_{1},w_{2})$ with $t_{2}$;

  3. 3.

read_adjacency($M,w,t$), which returns a list containing the values of $M$ in
the interval $[F(w,t,1),F(w,t,n)]$, i.e., the minimum possible timestamps to
arrive at any vertex while departing from $w$ at timestamp $t$;

  4. 4.

write_adjacency($M,w,t,L$), which replaces the values of $M$ values in the
interval $[F(w,t,1),F(w,t,n)]$ with the values of the list $L$, i.e., the
maximum possible timestamps to depart from any vertex while arriving at $w$ at
timestamp $t$.

Operations (1) and (2) access $O\left(1\right)$ pages on disk, while
operations (3) and (4) access $O\left(\nicefrac{{n}}{{B}}\right)$ pages, where
$B$ is the page size.

##  4 TTC Operations

n this section, we describe algorithms for the operations described in [4]:
the update operation add_contact(u,v,t); the query operations
can_reach($u$,$v$,$t_{1}$,$t_{2}$) and is_connected($t_{1}$,$t_{2}$); and the
reconstruction operation reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$). In
Section 4.1, we present our algorithm for add_contact(u,v,t) that receives a
contact and adds to our data structure the reachability information related to
the new available journeys passing thought it. In Section 4.2, we breafly
describe algorithms for can_reach($u$,$v$,$t_{1}$,$t_{2}$) and
is_connected($t_{1}$,$t_{2}$) since, as reachability information can be
directly accessed, they are straightforward. Finally, in Section 4.3, we
detail our algorithm for reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$) that
reconstructs a valid journey by concatenating one contact at a time.

###  4.1 Update operation

An algorithm to perform add_contact(u,v,t) must first add the reachability
information regarding the new trivial journey $\mathcal{J}_{triv}$ from vertex
$u$ to vertex $v$ departing at time $t$ and arriving at time $t+\delta$. Next,
for all vertices $w^{+}$ that $v$ can reach when departing at a time later
than or exactly $t+\delta$, the algorithm updates the reachability information
from $u$ to $w^{+}$ whether the new available journey passing through
$\mathcal{J}_{triv}$ has earlier arrival time. Then, for all vertices $w^{-}$
that can reach $u$ when arriving at a time earlier than or exactly $t$, the
algorithm updates the reachability information from $w^{-}$ to $v$ whether the
new available journey passing through $\mathcal{J}_{triv}$ has later departure
time. Finally, the algorithm must consider all new available journeys from
vertices $w^{-}$ to vertices $w^{+}$ that pass through $\mathcal{J}_{triv}$
and update the current reachability information if necessary.

Algorithm 1 describes the maintenance of both arrays $M_{out}$ and $M_{in}$
when inserting a new contact. In line 1, the algorithm checks if the structure
already has the information of the new contact $(u,v,t)$. If it still has not,
in line 2, it retrieves the latest departure timestamps of journeys departing
from vertices $w^{-}$ and arriving at vertex $u$ at timestamp $t$ as an array
$T^{-}$. In line 3, the algorithm retrieves the earliest arrival timestamps of
journeys departing from vertex $v$ at timestamp $t+\delta$ and arriving at
vertices $w^{+}$ as an array $T^{+}$. In lines 4 and 5, it sets the
reachability information about the new trivial journey
$\mathcal{J}_{triv}=(u,v,t)$ that departs at timestamp $t$ and arrives at
$t+\delta$. From line 6 to 14, the algorithm eagerly updates all cells
$M_{out}[w^{-},t^{\prime},w^{+}]=t^{+}$ for $t^{-}\geq t^{\prime}\geq 1$. In
this part, the algorithm proceeds by first iterating through all vertices
$w^{-}$, i.e., those that reached $u$ before than or exactly at timestamp $t$,
and retrieving their departure timestamps $t^{-}$. Then, it progressively
retrieves the current arrival timestamps to reach vertices $w^{+}$ when
departing at timestamp $t^{\prime}$, by reading the range
$[F_{out}(w^{-},t^{\prime},1),F_{out}(w^{-},t^{\prime},n)]$, and updates it
whether the new journeys passing through $\mathcal{J}_{triv}$ have earlier
arrival timestamps. Note that, vertices that could not reach $u$ before than
or exactly at timestamp $t$ have their arrival time equals to $-\infty$;
therefore, they are not considered in the while loop starting at line 8. This
process continues until the current reachability information in the whole
range does not change or $t^{\prime}<1$. Similarly, from line 6 to 14, the
algorithm eagerly updates all cells
$M_{in}[w^{+},t^{\prime\prime},w^{-}]=t^{-}$ for $t^{+}\leq
t^{\prime\prime}\leq\tau+\delta$. The algorithm proceeds by first iterating
through all vertices $w^{+}$, i.e., those that $v$ can reach departing after
or exactly at timestamp $t+\delta$, and retrieving their arrival timestamps
$t^{+}$. Then, it progressively retrieves the current departure timestamps in
which vertices $w^{-}$ departs present in range
$[F_{in}(w^{+},t^{\prime\prime},1),F_{in}(w^{+},t^{\prime\prime},n)]$ and
updates it whether the new available journeys passing through
$\mathcal{J}_{triv}$ have later departure timestamps. This process continues
until the current reachability information in the range does not change or
$t^{\prime\prime}>\tau+\delta$. Figure 2 illustrates the addition of new
contacts to a temporal graph along with the maintenance of the arrays
$M_{out}$ and $M_{in}$.

Algorithm 1 add_contact(u,v,t)

1: $u,v\in V$ with $u\neq{}v$, $n=|V|$, $t\in\mathcal{T}$, $\tau$, $\delta$,
$M_{out}$, $M_{in}$

2:if $\textsc{read\mathunderscore{}cell}(M_{out},u,t,v)\neq t+\delta$ then
$\triangleright$ check whether $(u,v,t)$ was inserted

3: $T^{-}\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{in},u,t)$

4:
$T^{+}\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{out},v,t+\delta)$

5: $T^{-}[u]\leftarrow t$ $\triangleright$ add the new trivial journey
information

6: $T^{+}[v]\leftarrow t+\delta$

7: for $w^{-}\text{ from }1\text{ up to }n$ do $\triangleright$ will update
$M_{out}$ with new journeys from $w^{-}$

8: $t^{\prime}\leftarrow T^{-}[w^{-}]$

9: while $t^{\prime}\neq-\infty$ and $t^{\prime}\geq 1$ do $\triangleright$
loop for $t^{-}\geq t^{\prime}\geq 1$

10:
$T^{+}_{cur}\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{out},w^{-},t^{\prime})$

11: $T^{+}_{cur}[w^{+}]\leftarrow\min(T^{+}_{cur}[w^{+}],T^{+}[w^{+}])$ for
$w^{+}\in[1,n]$

12: if $T^{+}_{cur}$ has not changed then

13: break

14:
$\textsc{write\mathunderscore{}adjacency}(M_{out},w^{-},t^{\prime},T^{+}_{cur})$

15: $t^{\prime}\leftarrow t^{\prime}-1$

16: for $w^{+}\text{ from }1\text{ up to }n$ do $\triangleright$ will update
$M_{in}$ with new journeys to $w^{+}$

17: $t^{\prime\prime}\leftarrow T^{+}[w^{+}]$

18: while $t^{\prime\prime}\neq\infty$ and $t^{\prime\prime}\leq\tau+\delta$
do $\triangleright$ loop for $t^{+}\leq t^{\prime\prime}\leq\tau+\delta$

19:
$T^{-}_{cur}\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{in},w^{+},t^{\prime\prime})$

20: $T^{-}_{cur}[w^{-}]\leftarrow\max(T^{-}_{cur}[w^{-}],T^{-}[w^{-}])$ for
$w^{-}\in[1,n]$

21: if $T^{-}_{cur}$ has not changed then

22: break

23:
$\textsc{write\mathunderscore{}adjacency}(M_{in},w^{+},t^{\prime\prime},T^{-}_{cur})$

24: $t^{\prime\prime}\leftarrow t^{\prime\prime}+1$

|  |  |
---|---|---|---
(a) added contact $(u,v,1)$

|  |  |
---|---|---|---
(b) added contacts $(u,w,3)$ and $(w,v,4)$

|  |  |
---|---|---|---
(c) added contact $(u,v,2)$

Figure 2:  Maintenance of our disk-based TTC, encoded as two arrays, $M_{out}$
(left table) and $M_{in}$ (right table), in different scenarios. Both arrays
are depicted as 2-dimensional arrays by grouping their first two dimensions.
In (a), the contact $(u,v,1)$ is inserted in the temporal graph and thus
$M_{out}$ and $M_{in}$ are updated using the information present in the left
and right expansions of the R-tuple $(u,v,1,2)$. In (b), both contacts
$(u,w,3)$ and $(w,v,4)$ are inserted and, additionally, a non-trivial journey
from $u$ to $v$ becomes possible. Finally, in (c), the contact $(u,v,2)$ is
inserted, allowing a faster journey departing from $u$ at time $2$ and
triggering the update of some cells of $M_{out}$ and $M_{in}$.

######  Theorem 2.

Algorithm 1 access $O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$ pages on disk.

###### Proof.

The read_cell operation in line 1 access a single page. The two read_adjacency
operations in lines 2 and 3 access $O\left(\nicefrac{{n}}{{B}}\right)$
sequential pages each. In lines 4 and 5, the algorithm writes the reachability
of the new trivial journey $J_{triv}=\\{(u,v,t)\\}$ in main memory. The for
loop starting at line 6 iterates over $n$ vertices $w^{-}$ and, the while loop
starting at line 8 iterates through $O(\tau)$ timestamps $t^{\prime}$. At each
of the $O(n\tau)$ iterations, it calls read_adjacency in order to read $n$
cells, and then (possibly) calls write_adjacency to write the $n$ cells back
while accessing, in each operation, $O\left(\nicefrac{{n}}{{B}}\right)$
sequential pages. Due to our mapping function $F_{out}$, at every timestamp
$t^{\prime}$, the algorithm will read a page that is arranged sequentially on
disk. The loop from line 15 to 23 does a similar computation. ∎

###  4.2 Reachability and Connectivity Queries

Both algorithms for can_reach($u$,$v$,$t_{1}$,$t_{2}$) and
is_connected($t_{1}$,$t_{2}$) are straightforward. The algorithm to perform
can_reach($u$,$v$,$t_{1}$,$t_{2}$) comprises testing whether
$\textsc{read\mathunderscore{}cell}(M_{out},u,t_{1},v)\leq{}t_{2}$ while
accessing only a single page from disk. The algorithm to perform
is_connected($t_{1}$,$t_{2}$), for each origin vertex $u\in V$, calls
$tmp\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{out},u,t_{1})$ and
then for each destination vertex $v\in V$, it checks whether
$tmp[v]\leq{}t_{2}$. As soon as a check is negative, the answer is false;
otherwise, it is true. Therefore, while the algorithm reads all cells, it
accesses $O\left(\nicefrac{{n^{2}}}{{B}}\right)$ pages on disk.

###  4.3 Journey Reconstruction

For the reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$) query, we need to augment
each cell of $M_{in}$ with the first successor vertex of the corresponding
journeys. Algorithm 1 can be trivially modified to include this information.
For instance, the successor vertex of a trivial journey from a contact
$(u,v,t)$ is the vertex $v$ since it is the first successor of $u$. Thus,
while composing previous R-tuples in our update algorithm, one would need to
read previous reachability information and compose than appropriately
considering also the successor vertex present in each cell of $M_{in}$.

Algorithm 2 gives the details to process the
reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$) query. Its goal algorithm is to
reconstruct a journey by unfolding the intervals and successor fields. In line
1, it initializes an empty journey $\mathcal{J}$. In line 2, it retrieves the
earliest timestamp $t^{+}$ a journey from vertex $u$ departing at timestamp
$t_{1}$ can arrive at vertex $v$ by reading on disk the entry
$M_{out}[u,t_{1},v]$. If $t^{+}\leq{}t_{2}$, it starts reconstructing the
resulting journey, otherwise, it returns an empty journey since there is no
journey completely in the interval $[t_{1},t_{2}]$. From lines 4 to 10, it
reconstructs the resulting journey by: first, in lines 4 and 5, initializing
the successor vertex $succ$ to $u$, and accessing on disk all the entries
$M_{in}[v,t^{+},w]$ for $w\in V$; then, from lines 6 to 10, the joursey s
reconstructed by iteratively accessing the next earliest departing timestamp
$t^{-}$ and the corresponding successor vertex $next\mathunderscore{}succ$
that reaches $v$ at timestamp $t^{+}$ while concatenating the contact
$(succ,next\mathunderscore{}succ,t^{-})$ at the end of $\mathcal{J}$ and
updating the current successor vertex.

Algorithm 2 reconstruct_journey($u$,$v$,$t_{1}$,$t_{2}$)

1:$[t_{1},t_{2}]\subset{\mathcal{T}},u,v\in V$ with $u\neq v$

2:$\mathcal{J}\leftarrow\\{\\}$

3:$(t^{+},\\_)\leftarrow\textsc{read\mathunderscore{}cell}(M_{out},u,t_{1},v)$

4:if $t^{+}\leq{}t_{2}$ then

5: $succ\leftarrow u$

6: $in\leftarrow\textsc{read\mathunderscore{}adjacency}(M_{in},v,t^{+})$

7: while $succ\neq v$ do

8: $t^{-}\leftarrow in[succ].t$

9: $next\mathunderscore{}succ\leftarrow in[succ].succ$

10:
$\mathcal{J}\leftarrow\mathcal{J}\cup(succ,next\mathunderscore{}succ,t^{-})$

11: $succ\leftarrow next\mathunderscore{}succ$

12:return $\mathcal{J}$

######  Theorem 3.

Algorithm 2 sequentially accesses $O\left(\nicefrac{{n}}{{B}}\right)$ pages on
disk, where $n$ is the number of vertices and $B$ is the page size.

###### Proof.

The algorithm accesses one page by calling
$\textsc{read\mathunderscore{}cell}(M_{out},u,t_{1},v)$ in line 2. After that,
it is known whether a journey exists or not. If a journey exists, it
sequentially accesses $\nicefrac{{n}}{{B}}$ pages by calling
$\textsc{read\mathunderscore{}adjacency}(M_{in},v,t^{+})$ in line 5. The
result of this call has all information needed to reconstruct a valid journey.
Finally, in the loop from line 6 to line 8, the algorithm extends the
resulting journey by one contact at each iteration using information already
in memory. Thus, the number of pages accessed is dominated by the call
$\textsc{read\mathunderscore{}adjacency}(M_{in},v,t^{+})$. ∎

##  5 Experiments

In this section, we will present experiments comparing our novel data
structure based on sequential arrays with the approach we adapted based on [4]
using B+-trees as a replacement for self-balanced binary search trees (BSTs).
Briefly, the approach introduced in [4] stores, in a matrix $n\times n$,
pointers to BSTs containing time intervals. In each BST, only non-redundant
intervals are kept, i.e., those that do not contain another interval in the
same tree. The authors proposed to use _join-based_ operations in order to
remove sequences of non-redundant intervals in $\log{\tau}$ time. These
operations can be found in LABEL:btreejoinsplit for B+-trees.

In the following, we present two experiments in Sections 5.1 and 5.2. In the
first one, we inserted unsorted contacts from complete temporal graphs,
incrementally, in both data structures using the operation add_contact(u,v,t).
In the second one, we inserted shuffled contacts from real-world datasets.

###  5.1 Experiments with Synthetic Data

In this first experiment, we generated complete temporal graphs with the
number of vertices fixed to $100$ and varied the number of timestamps $\tau$
from $10$ to $10000$. Then, we inserted their shuffled contacts in both data
structures using the add_contact(u,v,t) operation. The time to preallocate and
initialize the arrays on disk for our array-based data structure was not
considered in the total time. We note that this extra cost can be high for
large parameters; therefore, one should consider it whenever applicable.

Figure 3: Cumulative wall-clock time to maintain data structures for
reachability queries on synthetic data. We inserted shuffled contacts from
complete temporal graphs into the data structures varying the number of
timestamps $\tau$ while fixing the number of vertices to $100$. Red lines
represent our novel data structure based on sequential arrays. Blue lines
represent our adaptation of the approach introduced in [4] using B+-trees as
self-balanced BSTs.

Figure 3 shows the mean cumulative wall-clock time, averaged over $10$
executions, to maintain both data structures as new unsorted contacts were
inserted. We see that our novel data structure performs better for all
configurations. Even though the worst-case complexity of our algorithm for the
add_contact(u,v,t) operation is linear in $\tau$ instead of logarithmic, it
runs much faster using synthetic data. We attribute this behavior to the fact
that as new contacts are inserted, the probability of composing better
R-tuples, i.e. journeys, decreases rapidly and, thus, our data structure
updates on average only few cells per contact insertion.

Next, we will argue why the run time of our algorithm reduces with the
addition of contacts. Each pair of vertices $(u,v)$ are associated to a set
$\mathcal{I}$ containing intervals $[t^{-},t^{+}]\subseteq[1,\tau]$ in which
$u$ can reach $v$ departing at $t^{-}$ and arriving at $t^{+}$. For a
particular pair of vertices, when an algorithm inserts a new interval $I$, all
intervals $I^{\prime}$ such that $I\subseteq I^{\prime}$ can be safely removed
since they become redundant. Our data structure organizes these intervals in
the arrays $M_{out}$ and $M_{in}$, which fix, respectively, the left and right
endpoints, and our update algorithm discards redundant intervals by updating
their cells accordingly by using Definition 6.

Consider the hierarchy of intervals illustrated in Figure 4(a) for $\tau=4$.
Each interval with length $l$ is linked to the intervals with length $l-1$
that it totally encloses. For example, interval $[0,5]$, with length $5$,
links to intervals $[0,4]$ and $[1,5]$, with length $4$, because
$[0,4]\subseteq[0,5]$ and $[1,5]\subseteq[0,5]$. Initially, all intervals are
available for insertion in our data structure. When a new interval $[1,2]$ is
inserted, as show in Figure 4(b), all intervals that contain it, including
itself, are not available for insertion anymore.

Our update algorithm conceptually removes these intervals by drawing left and
right frontiers separating available and non-available intervals starting from
$[1,2]$. For instance, intervals $[1,2]$ and $[0,2]$, which belong to the left
frontier, are updated in $M_{in}$ since they share the same right endpoint,
and intervals $[1,2]$, $[1,3]$, $[1,4]$ and $[1,5]$, which belong to the right
frontier, are updated in $M_{out}$ since they share the same left endpoint. In
this proccess, up to $\tau$ cells are updated in both $M_{in}$ and $M_{out}$.

Next, when a new interval $[3,5]$ is inserted, as shown in Figure 4(c), our
algorithm must, again, draw the left and right frontiers starting from
$[3,5]$; however it does not need to advance previously drawn frontiers. In
this case, only intervals $[3,5]$ and $[2,5]$, which belong to the left
frontier, are updated in $M_{in}$. In Figure 4(d), interval $[2,3]$ is
inserted and the same proccess repeats. We see that as new intervals are
inserted, the number of available intervals rapidly reduces. Thus, even though
our algorithm has complexity $O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$, it
can run much faster when considering a sequence of contact insertions since
the number of cells to be updated reduces over time.

Moveover, it is guaranteed that, for each new contact $(u,v,t)$, our algorithm
will make unavailable for insertion every interval that is still available
inside and at the frontiers starting from $[t,t+\delta]$ in the lowest level
of the hierarchy associated with $(u,v)$.

|
---|---
(a) all available intervals | (b) insert $[1,2]$
|
(c) insert $[3,5]$ | (d) insert $[2,3]$
Figure 4: Illustration of the proccess performed by our update algorithm
considering a fixed pair of vertices $(u,v)$ from a temporal graph with
$\tau=4$. Available intervals for insertion are colored in black, and
invalidated intervals, i.e. intervals that should not be considered anymore by
our update algorithm, are colored in different colors. Links represent the
direct containement relation between intervals with length $l$ and intervals
with length $l-1$.

###  5.2 Experiments with Real-World Datasets

In this second experiment, we downloaded small and medium real-world available
on https://networkrepository.com/dynamic.php, and preprocessed them using our
script available on https://bitbucket.com/luizufu/temporalgraph-datasets-
preprocessing. During the preprocessing, we rellabelled the vertices and
shifted the timestamps of each dataset so that vertex identifiers were
beetween $[1,n]$ and timestamp values start from $1$. Then, we inserted the
shuffled contacts of each dataset in both data structures using the
add_contact(u,v,t) operation. We assumed that all used datasets represent
temporal digraphs, and we used $\delta=1$, i.e., traversing any contact takes
one time unit.

dataset | $n$ | $\tau$ | contacts | density | Array-Based | Tree-Based
---|---|---|---|---|---|---
aves-sparrow | $52$ | $2$ | $516$ | $0.1$ | $0.01\pm 0$ | $0.07\pm 0$
aves-weaver | $445$ | $23$ | $1423$ | $0.003$ | $0.19\pm 0$ | $1.16\pm 0.01$
aves-wildbird | $202$ | $6$ | $11900$ | $0.05$ | $0.97\pm 0.01$ | $9.52\pm 0.15$
ant-colony1 | $113$ | $41$ | $111578$ | $0.46$ | $25.84\pm 0.19$ | $161.3\pm 1.03$
ant-colony2 | $131$ | $41$ | $139925$ | $0.2$ | $43.96\pm 0.49$ | $261.98\pm 1.96$
ant-colony3 | $160$ | $41$ | $241280$ | $0.23$ | $89.37\pm 0.73$ | $524.79\pm 6.71$
ant-colony4 | $102$ | $41$ | $81599$ | $0.19$ | $16.49\pm 0.12$ | $104.62\pm 1.27$
ant-colony5 | $152$ | $41$ | $194317$ | $0.21$ | $166.93\pm 3.63$ | $526.03\pm 144.32$
ant-colony6 | $164$ | $39$ | $247214$ | $0.24$ | $88.99\pm 0.93$ | $608.87\pm 167.97$
copresence-LH10 | $73$ | $259181$ | $150126$ | $0.0001$ | - | $61.5\pm 0.53$
copresence-LyonSchool | $242$ | $117721$ | $6594492$ | $0.001$ | - | $14887.45\pm 1576.49$
kilifi-within-households | $54$ | $59$ | $32426$ | $0.19$ | $0.09\pm 0$ | $0.21\pm 0$
mammalia-primate | $25$ | $19$ | $1340$ | $0.12$ | $0.09\pm 0$ | $0.33\pm 0.01$
mammalia-raccoon | $24$ | $52$ | $1997$ | $0.06$ | $0.21\pm 0$ | $0.44\pm 0.01$
mammalia-voles-bhp | $1686$ | $63$ | $5324$ | $0.00003$ | $13.76\pm 0.82$ | $19.22\pm 0.24$
mammalia-voles-kcs | $1218$ | $64$ | $4258$ | $0.00004$ | $7.92\pm 0.27$ | $11.78\pm 0.09$
mammalia-voles-plj | $1263$ | $64$ | $3863$ | $0.00003$ | $6.18\pm 0.23$ | $10.68\pm 0.04$
mammalia-voles-rob | $1480$ | $63$ | $4569$ | $0.00003$ | $10.28\pm 0.41$ | $15.09\pm 0.12$
tortoise-bsv | $136$ | $4$ | $554$ | $0.008$ | $0.01\pm 0$ | $0.14\pm 0.01$
tortoise-cs | $73$ | $10$ | $258$ | $0.005$ | $0.01\pm 0$ | $0.05\pm 0$
tortoise-fi | $787$ | $9$ | $1713$ | $0.0003$ | $0.15\pm 0$ | $2.71\pm 0.01$
trophallaxis-colony1 | $41$ | $8$ | $308$ | $0.02$ | $0.02\pm 0$ | $0.06\pm 0$
trophallaxis-colony2 | $39$ | $8$ | $330$ | $0.03$ | $0.02\pm 0$ | $0.05\pm 0$
Table 2: Total wall-clock time in seconds to insert all shuffled contacts from
real-world datasets with number of vertices $n$, number of timestamps $\tau$,
number of contacts into data structures for reachability queries, and the
density of the temporal graph represented by the dataset. Values were rounded
to two decimal places. Array-based refers to our novel data structure and
tree-based refers to our implementation of the approach introduced in [4]
using B+-trees as BSTs replacement. Executions that reached the time limit of
$5$ hours are marked with the symbol “-”.

Table 3 shows the mean wall-clock time, averaged over $10$ executions, to
insert all shuffled contacts of each dataset into both data structures. We see
that our novel data structure performs better on the majority of datasets.
However for the largest datasets, copresence-LH10 and copresence-LyonSchool,
the tree-based data structure performed better. Both datasets have high values
for $\tau$ and low density. It means that, as density is too small, each
insertion of a contact $(u,v,t)$ may trigger an initial update over arrays
$M_{out}$ and $M_{in}$ that will touch many cells on disk. As in Figure 4(b),
for most insertions, our update algorithm will draw left and right frontiers
on the almost empty hierarchy associated with the pair of vertices $(u,v)$
starting from interval $[t,t+\delta]$. Therefore, in this case, the linear
factor on $\tau$ from the cost $O(\nicefrac{{n^{2}\tau}}{{B}})$ of our update
algorithm will have a bigger impact on the run time since the sequence of
insertions is not sufficiently long for our algorithm to benefit from later
insertions.

##  6 Concluding remarks

We presented in this paper an incremental disk-based data structure to solve
the dynamic connectivity problem in temporal graphs. Our data structure
prioritizes query time, answering reachability queries by accessing only one
page. Based on the ability to quickly retrieve reachability information among
vertices inside time intervals, it can: insert contacts in a non-chronological
order accessing $O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$ pages, where $B$
is the size of disk pages; check whether a temporal graph is connected within
a time interval accessing $O\left(\nicefrac{{n^{2}}}{{B}}\right)$ pages, and
reconstruct journeys accessing $O\left(\nicefrac{{n}}{{B}}\right)$ pages. Our
algorithms exploit the special features of non-redundant (minimal)
reachability information, which we represent explicitly through the concept of
expanded R-tuples. As in [4], the core of our data structure, is essentially a
collection of non-redundant R-tuples, whose size (and that of the data
structure itself) cannot exceed $O\left(n^{2}\tau\right)$. However, in our
approach, all this space must be preallocated on disk. The benefit of our data
structure is that algorithms explicitly manage data sequentially and,
therefore, it is more suitable for secondary memories in which random accesses
are expensive.

Further investigations could be done in the direction of improving the
complexity of our update algorithm. Can add_contact(u,v,t) access less than
$O\left(\nicefrac{{n^{2}\tau}}{{B}}\right)$ pages? Another direction could be
designing efficient disk-based data structures for the decremental and the
fully-dynamic versions of this problem. With _unsorted_ contact insertion and
deletion, it seems to represent both a significant challenge and a natural
extension of the present work, one that would certainly develop further our
common understanding of temporal reachability. Finally, it could be worth to
investigate compressing algorithms to reduce the space of our data structure
and the number of pages accessed by our update algorithm. Specifically, we
think that compression algorithms based on differences and run-length coding
[10] could achieve a very high compression rate since the arrays $M_{out}$ and
$M_{in}$ store repeating ordered values. The compressing schema could also
solve the preallocation and initialization problem since all cells of
$M_{out}$ and $M_{in}$ have, initially, the same value, which are very
compressible.

##  Appendix A Join and split operations for B+-trees

Let each leaf node of B+-trees contains an array $K$ of keys of size $N$ and a
pointer to its next sibling. Let each non-leaf node contains an array of keys
$K$ of size $M$ and an additional array of pointers $C$ to child nodes of size
$M+1$. Additionally, assume that every node contains the height of the sub-
tree it belongs, a pointer to its leftmost leaf child and a pointer to its
rightmost leaf child. We note that we use these additional per-node data to
simplify our algorithms and discussions. In a real implementation, only the
root node (the tree itself) must maintain them during the insertion and update
operations. Information regarding the rest of the nodes can be computed during
the execution of the next algorithms without increasing complexities.

###  A.1 Join operation on B+-trees

Algorithm 3 performs the operation join for B+-trees. Given two B+-trees
$T_{left}$ and $T_{right}$, such that keys present in $T_{left}$ are smaller
to keys in $T_{right}$, it must merge both trees in order to create a new
valid B+-tree $T$ containing all keys present in $T_{left}$ and $T_{right}$.
As B+-trees place leaf nodes at the same level, it simply inserts or shares
the data present in the root node of the smaller tree into the appropriate
node at the same height in the bigger tree. Then it maintains the B+-tree
invariances up to its root node of the changed bigger tree whenever necessary.
First, in line 11, the algorithm sets the next sibling of the rightmost leaf
of $T_{left}$ to be the leftmost leaf of $T_{right}$. Then, if
$\texttt{height}(T_{left})\geq\texttt{height}(T_{right})$, in lines 13 and 14,
it adds $T_{right}$ to $T_{left}$ by calling
$\textsc{joinRight}(T_{left},T_{right})$ and returns $T_{left}$; otherwise, in
lines 16 and 17, it adds $T_{left}$ to $T_{right}$ by calling
$\textsc{joinLeft}(T_{left},T_{right})$ and returns $T_{right}$. From lines 1
to 11, we detail the joinRight routine, the joinLeft routine is implemented
symmetrically. In line 1, the algorithm descends $T_{left}$ until reaching the
rightmost node $n_{left}$ at the same height of the root node of $T_{right}$.
If a single node of size $B$ can fit the content of both $n_{left}$ and the
root node of $T_{right}$, in line 5, it simply merges both nodes by adding to
$n_{left}$ the data present in the root node of $T_{right}$. Otherwise, in
line 7, it equally shares the data of both nodes, and, in line 8, it inserts
into the parent of $n_{lelft}$ a new key together with a pointer to the root
node of $T_{right}$. If the parent node has no space left to accommodate the
new data, a node splitting routine must be invoked and this process can
continue up to the root node of $T_{left}$. Finally, if the algorithm needs to
split the current root node of $T_{left}$, then it creates a new root node,
and, in this case, In line 10, it increments the height of $T_{left}$ by one.

Algorithm 3 join

1:Two trees of intervals $T_{left}$ and $T_{right}$

2:procedure joinRight($T_{left},T_{right}$)

3:
$n_{left}\leftarrow\texttt{descend\mathunderscore{}right}(T_{left},\texttt{height}(T_{left})-\texttt{height}(T_{right}))$

4: $n_{right}\leftarrow\texttt{root}(T_{right})$

5: if $\texttt{size}(n_{left})+\texttt{size}(n_{right})\leq B$ then

6: $n_{left}\leftarrow\texttt{merge}(n_{left},n_{right})$

7: else

8: $\texttt{share}(n_{left},n_{right})$

9:
$\texttt{insert\mathunderscore{}rec}(\texttt{parent}(n_{left}),\texttt{min\mathunderscore{}key}(n_{right}),n_{right})$

10: if new root node was created then

11: $\texttt{height}(T_{left})\leftarrow\texttt{height}(T_{left})+1$

12:$\texttt{next\mathunderscore{}leaf}(\texttt{rightmost\mathunderscore{}leaf}(T_{left}))\leftarrow\texttt{leftmost\mathunderscore{}leaf}(T_{right})$

13:if $\texttt{height}(T_{left})\geq\texttt{height}(T_{right})$ then

14: $\textsc{joinRight}(T_{left},T_{right})$

15: return $T_{left}$

16:else

17: $\textsc{joinLeft}(T_{left},T_{right})$

18: return $T_{right}$

######  Theorem 4.

Algorithm 3 accesses
$O(|\texttt{height}(T_{left})-\texttt{height}(T_{right})|)$ pages in the
worst-case.

###### Proof.

In line 11, the algorithm accesses one page to set the next child node of the
rightmost leaf node of $T_{left}$. Then, it calls joinRight or joinLeft
depending on the heights of $T_{left}$ and $T_{right}$, both accessing the
same amount of pages. Without loss of generality, assume that
$\texttt{height}(T_{left})\geq\texttt{height}(T_{right})$ and it calls thus
the joinRight procedure. Then, at line 2, the algorithm accesses
$\texttt{height}(T_{left})-\texttt{height}(T_{right})$ pages while descending
to the rightmost node of $T_{left}$ at height $\texttt{height}(T_{right})$.
Next, if there is enough room to fit the data of both nodes being merged in a
single node, it accesses $O(1)$ pages and the algorithm ends. Otherwise, it
accesses $O(1)$ pages to share the content present in the considered nodes.
Then, it accesses, again,
$O(\texttt{height}(T_{left})-\texttt{height}(T_{right}))$ pages in order to
insert new key and pointer pairs up to the root of $T_{left}$ in the worst-
case. Finally, if a new node is created, it access one more page to increment
the height of $T_{left}$ and the algorithm ends. ∎∎

###  A.2 Split operation on B+-trees

Algorithm 4 performs the operation split for B+-trees. Given an interval key
$L$, it must split a tree $T$ in two trees $T_{left}$ and $T_{right}$ such
that all keys in $T_{left}$ are smaller than $I$ and all keys in $T_{right}$
are greater or equal to $I$. To accomplish this task, it recursively descends
$T$ from the root node to the leaf node containing the biggest key less than
$L$ while partitioning nodes appropriately and, during the backward phase of
the recursion, progressively building $T_{left}$ and $T_{right}$. During each
recursive step, in line 1, the algorithm first finds the position $k$ in the
current root node such that $K[k]\geq I$, where $C[k]$ is the pointer that
branches to the next child node for non-leaf nodes. If the current root node
is a leaf, in line 3, it partitions the current node in two sub-trees:
$T_{left}$, containing a node with $K[1\ldots k-1]$; and $T_{right}$,
containing a node with $K[k\ldots N]$. Then, in line 4, it sets the next
sibling of $T_{left}$’s root to $nil$; in line 5, it sets the next sibling of
$T_{right}$ to the next sibling of $T$’s root; and, in line 6, it returns
$(T_{left},T_{right})$. Note that no other leaf node besides the affected ones
must update the pointer to its next siblings since the resulting trees will
reuse the previous linkages. Next, if the current node is a non-leaf, in line
7, the algorithm partitions the current node in three sub-trees: $T_{left}$,
containing a root node with $K[1\ldots k-2]$ and $C[1\ldots k-1]$;
$T_{child}$, containing a root node with $K[k-1\ldots k]$ and $C[k]$; and (3)
$T_{right}$, containing a root node with $K[k+1\ldots M]$ and $C[k+1\ldots
M+1]$. Additionally, whenever a sub-tree have only one pointer in its root
node, its respective root node becomes the child pointed by it in order to
maintain the correct B+-tree layout. Then, in line 8, it calls the split
algorithm itself passing $T_{child}$ as parameter and obtaining two sub-trees
$T^{\prime}_{left}$ and $T^{\prime}_{right}$ as the intermediate result.
Finally, in line 9, it advances the intermediate result by joining them
appropriately with the sub-trees of the current level.

Algorithm 4 split

1:A tree of intervals $T$ and a key interval L

2:$k\leftarrow\texttt{find\mathunderscore{}key\mathunderscore{}position}(\texttt{root}(T),L)$

3:if $\texttt{root}(T)$ is a leaf then

4:
$(T_{left},T_{right})\leftarrow\texttt{split\mathunderscore{}leaf}(root(T),k)$

5: $\texttt{next\mathunderscore{}leaf}(\texttt{root}(T_{left}))\leftarrow nil$

6:
$\texttt{next\mathunderscore{}leaf}(\texttt{root}(T_{right}))\leftarrow\texttt{next\mathunderscore{}leaf}(\texttt{root}(T))$

7: return $(T_{left},T_{right})$

8:$(T_{left},T_{child},T_{right})\leftarrow\texttt{split\mathunderscore{}non\mathunderscore{}leaf}(root(T),k)$

9:$(T^{\prime}_{left},T^{\prime}_{right})\leftarrow\texttt{split}(T_{child},L)$

10:return
$(\texttt{join}(T_{left},T^{\prime}_{left}),\texttt{join}(T^{\prime}_{right},T_{right}))$

######  Theorem 5.

Algorithm 4 accesses $O(\log_{B}{(\tau)})$ pages in the worst-case where
$\tau$ is the maximum number of keys in the tree.

###### Proof.

During each recursive step, Algorithm 4 needs to: (1) partition the current
node being considered in at least two sub-trees; (2) change pointers to next
leaf siblings, whether the current node is a leaf; and (3) join the current
sub-trees with the sub-trees resulting from the next recursive step. For (1),
the algorithm accesses $O(1)$ pages. In the worst-case scenario, if the root
node of a sub-tree has only a single child, the algorithm reads this child in
order to make it the new root node. For (2), the algorithm also accesses
$O(1)$ pages since only two leaf nodes are updated. For (3), the algorithm
calls the join algorithm twice. Without loss of generality, consider only
calls maintaining $T_{left}$. From the node above the leaf node containing the
split key up to the root of $T$, the algorithm joins the current left sub-tree
$T_{left}$ with the intermediate left sub-tree $T^{\prime}_{left}$ resulting
from the previous iteration. At each iteration there is a
$\texttt{join}(T_{left},T^{\prime}_{left})$ call in which $T_{left}$ is either
empty or non-empty. In case $T_{left}$ is empty, the algorithm pays nothing
and, at the next iteration, the difference in height between $T_{left}$ and
$T^{\prime}_{left}$ increases by one. In case $T_{left}$ is non-empty, the
algorithm pays the difference in height accumulated so far and, at the next
iteration, the difference in height resets to one. Therefore, as the summation
of all payments is at most the height of the tree, the algorithm accesses
$O(\log_{B}(\tau))$ pages while processing all join calls. ∎∎

#### Acknowledgements

This study was financed in part by Fundação de Amparo à Pesquisa do Estado de
Minas Gerais (FAPEMIG) and the Coordenação de Aperfeiçoamento de Pessoal de
Nível Superior - Brasil (CAPES) - Finance Code 001* - under the “CAPES PrInt
program” awarded to the Computer Science Post-graduate Program of the Federal
University of Uberlândia.

## References

  * [1] David J Abel  “A B+-tree structure for large quadtrees”  In _Computer Vision, Graphics, and Image Processing_ 27.1 Elsevier, 1984, pp. 19–31
  * [2] Matthieu Barjon et al.  “Testing temporal connectivity in sparse dynamic graphs”  In _CoRR_ abs/1404.7634, 2014
  * [3] L. Bedogni, M. Fiore and C. Glacet  “Temporal reachability in vehicular networks”  In _IEEE INFOCOM 2018 - IEEE Conference on Computer Communications_ , 2018, pp. 81–89
  * [4] Luiz F. Brito, Marcelo Albertini, Arnaud Casteigts and Bruno A.. Travençolo  “A Dynamic Data Structure for Temporal Reachability with Unsorted Contact Insertions”, 2021  arXiv:2102.04187 [cs.DS]
  * [5] Daniel Bryce and Subbarao Kambhampati  “A tutorial on planning graph based reachability heuristics”  In _AI Magazine_ 28.1, 2007, pp. 47
  * [6] Leo Cacciari and Omar Rafiq  “A temporal reachability analysis”  In _Protocol Specification, Testing and Verification XV: Proceedings of the Fifteenth IFIP WG6.1 International Symposium on Protocol Specification, Testing and Verification, Warsaw, Poland, June 1995_ Boston, MA: Springer US, 1996, pp. 35–49
  * [7] Arnaud Casteigts, Paola Flocchini, Walter Quattrociocchi and Nicola Santoro  “Time-varying graphs and dynamic networks”  In _International Journal of Parallel, Emergent and Distributed Systems_ 27.5 Taylor & Francis, 2012, pp. 387–408
  * [8] Andrea E.F. Clementi et al.  “Flooding Time in Edge-Markovian Dynamic Graphs”  In _Proceedings of the Twenty-Seventh ACM Symposium on Principles of Distributed Computing_ , PODC ’08  Toronto, Canada: Association for Computing Machinery, 2008, pp. 213–222
  * [9] Edith Cohen, Eran Halperin, Haim Kaplan and Uri Zwick  “Reachability and Distance Queries via 2-Hop Labels”  In _SIAM Journal on Computing_ 32.5, 2003, pp. 1338–1355
  * [10] Patrick Damme, Dirk Habich, Juliana Hildebrandt and Wolfgang Lehner  “Lightweight Data Compression Algorithms: An Experimental Survey ( Experiments and Analyses)”  In _EDBT_ , 2017
  * [11] Betsy George, Sangho Kim and Shashi Shekhar  “Spatio-temporal network databases and routing algorithms: a summary of results”  In _Advances in Spatial and Temporal Databases_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 460–477
  * [12] Haixun Wang et al.  “Dual Labeling: Answering Graph Reachability Queries in Constant Time”  In _22nd International Conference on Data Engineering (ICDE’06)_ , 2006, pp. 75–75
  * [13] Khandaker Tabin Hasan, Sheak Rashed Haider Noori, Abdus Salam and Md Anwarul Kabir  “Making sense of time: timeline visualization for public transport schedule”  In _Symposium on Human-Computer Interaction and Information Retrieval (HCIR 2011), Washington_ , 2011
  * [14] C. Hurter et al.  “Bundled visualization of dynamicgraph and trail data”  In _IEEE Transactions on Visualization and Computer Graphics_ 20.8, 2014, pp. 1141–1157
  * [15] Alexandre Camargo Martensen, Santiago Saura and Marie-Josee Fortin  “Spatio-temporal connectivity: assessing the amount of reachable habitat in dynamic landscapes”  In _Methods in Ecology and Evolution_ 8.10, 2017, pp. 1253–1264
  * [16] Vera Zaychik Moffitt and Julia Stoyanovich  “Portal: A Query Language for Evolving Graphs”  In _CoRR_ abs/1602.00773, 2016  arXiv: http://arxiv.org/abs/1602.00773
  * [17] Sebastiaan J. Schaik and Oege Moor  “A Memory Efficient Reachability Data Structure through Bit Vector Compression”  In _Proceedings of the 2011 ACM SIGMOD International Conference on Management of Data_ , SIGMOD ’11  Athens, Greece: Association for Computing Machinery, 2011, pp. 913–924
  * [18] S. Seufert, A. Anand, S. Bedathur and G. Weikum  “FERRARI: Flexible and efficient reachability range assignment for graph indexing”  In _2013 IEEE 29th International Conference on Data Engineering (ICDE )_ , 2013, pp. 1009–1020
  * [19] John Tang, Mirco Musolesi, Cecilia Mascolo and Vito Latora  “Characterising temporal distance and reachability in mobile and online social networks”  In _SIGCOMM Comput. Commun. Rev._ 40.1 New York, NY, USA: Association for Computing Machinery, 2010, pp. 118–124
  * [20] Hao Wei, Jeffrey Xu Yu, Can Lu and Ruoming Jin  “Reachability querying: an independent permutation labeling approach”  In _The VLDB Journal_ 27.1 Springer, 2018, pp. 1–26
  * [21] John Whitbeck, Marcelo Dias de Amorim, Vania Conan and Jean-Loup Guillaume  “Temporal reachability graphs”  In _Proceedings of the 18th Annual International Conference on Mobile Computing and Networking_ , Mobicom ’12  Istanbul, Turkey: Association for Computing Machinery, 2012, pp. 377–388
  * [22] Matthew J. Williams and Mirco Musolesi  “Spatio-temporal networks: reachability, centrality and robustness”  In _Royal Society Open Science_ 3.6, 2016, pp. 160196
  * [23] G. Wu et al.  “Mining spatio-temporal reachable regions over massive trajectory data”  In _2017 IEEE 33rd International Conference on Data Engineering (ICDE )_ , 2017, pp. 1283–1294
  * [24] H. Wu et al.  “Reachability and time-based path queries in temporal graphs”  In _2016 IEEE 32nd International Conference on Data Engineering (ICDE )_ , 2016, pp. 145–156
  * [25] W. Zeng et al.  “Visualizing mobility of public transportation system”  In _IEEE Transactions on Visualization and Computer Graphics_ 20.12, 2014, pp. 1833–1842
  * [26] Andy Diwen Zhu, Wenqing Lin, Sibo Wang and Xiaokui Xiao  “Reachability Queries on Large Dynamic Graphs: A Total Order Approach”  In _Proceedings of the 2014 ACM SIGMOD International Conference on Management of Data_ , SIGMOD ’14  Snowbird, Utah, USA: Association for Computing Machinery, 2014, pp. 1323–1334