# Type III solar radio burst detection and classification: A deep learning
approach

Jeremiah Scully  Dept. of Computer and Software Engineering
Athlone Institute of Technology
Athlone, Ireland
<EMAIL_ADDRESS>Ronan Flynn  Dept. of Computer and Software
Engineering
Athlone Institute of Technology
Athlone, Ireland
<EMAIL_ADDRESS>Eoin Carley  School of Cosmic Physics
Dublin Institute of Advanced Studies
Dublin, Ireland
<EMAIL_ADDRESS>Peter Gallagher  School of Cosmic Physics
Dublin Institute of Advanced Studies
Dublin, Ireland
<EMAIL_ADDRESS>Mark Daly  Dept. of Computer and Software Engineering
Athlone Institute of Technology
Athlone, Ireland
<EMAIL_ADDRESS>
###### Abstract

Solar Radio Bursts (SRBs) are generally observed in dynamic spectra and have
five major spectral classes, labelled Type I to Type V depending on their
shape and extent in frequency and time. Due to their complex characterization,
a challenge in solar radio physics is the automatic detection and
classification of such radio bursts. Classification of SRBs has become
fundamental in recent years due to large data rates generated by advanced
radio telescopes such as the LOw-Frequency ARray, (LOFAR). Current state-of-
the-art algorithms implement the Hough or Radon transform as a means of
detecting predefined parametric shapes in images. These algorithms achieve up
to 84% accuracy, depending on the Type of radio burst being classified. Other
techniques include procedures that rely on Constant-False-Alarm-Rate
detection, which is essentially detection of radio bursts using a de-noising
and adaptive threshold in dynamic spectra. It works well for a variety of
different Types of radio bursts and achieves an accuracy of up to 70%. In this
research, we are introducing a methodology named You Only Look Once v2
(YOLOv2) for solar radio burst classification. By using Type III simulation
methods we can train the algorithm to classify real Type III solar radio
bursts in real-time at an accuracy of 82.63% with a maximum 77 frames per
second (fps).

##  I Introduction

Ever since astronomer Richard Christopher Carrington observed the first ever
solar flare back in 1859, researchers have been both astonished and perplexed
by the behavior of the suns weather and in particular solar flares.

Solar flares are the most energetic explosive phenomena in the solar system,
often involving the acceleration of particles to near light-speed[1]. The
accelerated particles emit light from across the entire electromagnetic
spectrum, from gamma rays to radio waves. Radio emission is often high-
intensity and is generally observed as complex patterns in dynamic spectra
known as solar radio bursts[2]. SRBs generally come in five classifications,
labelled Type I to Type V, depending on their shape in dynamic spectra. They
can occur at a rate of thousands per day (especially Type III bursts), so it
is a computational challenge to automatically detect their occurrence and
determine their spectral characteristics. However, due to their complex
characterisation, classifying these radio bursts is a fundamental challenge.
This challenge has been made more complex in recent years with new technology
such as LOFAR providing high-volume data streams (up to 3 Gb/s at a single
station) of radio burst observations that need to be classified with high
speed and accuracy. The necessity of automated data pipelines for solar radio
bursts has been made even more apparent in recent years, with the the design
of LOFAR for Space Weather (LOFAR4SW)[3], a system upgrade which aims at
autonomously monitoring of solar radio activity (as well as heliospheric and
ionospheric activity). Software pipelines to automatically detect SRBs will be
an indispensable part of such a system in the near future. This paper
establishes the role of machine learning playing a crucial part of such a
pipeline.

There have been various attempts to automatically detect SRBs in dynamic
spectra in recent years. The current state-of-the-art involve algorithms that
implement the Hough or Radon transform as a means of detecting predefined
parametric shapes in images[4]. These algorithms achieve up to 84% accuracy,
depending on the Type of radio burst being classified. Other techniques
include procedures that rely on Constant-False-Alarm-Rate (CRAF) detection[5],
which is essentially detection of radio bursts using a de-noising and adaptive
threshold in dynamic spectra. It works well for a variety of different Types
of radio bursts and achieves an accuracy of up to 70%.

In recent years neural networks have been applied to the problem, in which
multi-modal deep learning was applied to a spectrogram at millimetric
wavelengths[6]. The system used auto-encoders and standard regularization in
tandem to achieve a burst detection accuracy of 82%, however this has not been
applied to metric wavelengths (the range of LOFAR) where the bursts can have
much more complex shapes.

Recent research has turned to object detection algorithms such as Faster R-CNN
to identify solar radio bursts[7]. This deep learning neural network proved to
be accurate at extracting small features of solar radio bursts with an average
precision (AP) of 91% however, it doesn’t have the performance of real time
detection. One other area in which machine learning is working in tandem with
radio telescope observing is SETI (Search for Extra-Terrestrial
Intelligence)[8]. SETI use their Allen telescope array to observe planetary
systems searching for Fast Radio Bursts (FRBs) using machine learning. SETI
highlights regular noise frequency’s and Radio Frequency Interference (RFI)
and then isolates irregular high frequency FRB spikes within the dynamic
spectrum using a deep Convolutional Neural Network (CNN) called a ResNet[9].
This model produced a recall score of 95%.

With recent research turning to deep CNNs and object detection for classifying
and detecting radio frequencies we decided to explore these areas even
further. There are many different forms of CNNs for object detection such as
You Only Look Once (YOLO)[10], Single Shot Detectors[11], Region-CNN
(R-CNN)[12], Fast R-CNN [13], Faster R-CNN[14] and Mask R-CNN[15]. Although
the methods mentioned have been proven to be very successful for object
detection YOLO has been the only algorithm to offer high accuracy and real
time detections on some datasets.

In this research, we test the accuracy of the deep learning algorithm
YOLOv2[16]. using the darkflow framework when applied to Type III solar radio
bursts. Using methods to simulate Type III bursts we can create a training set
of 80,000 images and train YOLOv2 to identify Type III radio bursts within a
dynamic spectrum at an accuracy of 82.63% at real-time frame rates of above
60fps.
In section II of this paper, we will be discussing what LOFAR is and what SRBs
are along with the dataset and model configuration used for this research. In
section III, we will discuss the training and test set, and also YOLOv2s
performance and accuracy on Type III data.

##  II Methodology

###  II-A LOFAR and I-LOFAR

LOFAR is a radio interferometer constructed in the north of the Netherlands
and across Europe which includes Ireland’s very own I-LOFAR station. LOFAR
covers the largely unexplored low frequency range from 10–240 MHz and offers a
unique number of observing capabilities. LOFAR makes use of Digital beam-
forming techniques meaning LOFAR can rapidly re-point it’s view to other
points of interest in a short space of time. It can also operate
simultaneously using multiple stations. LOFAR station can work both in tandem
with each other as one whole large radio telescope or singularly as a
standalone station. LOFAR antenna stations provide the exact same simple
functions as radio dishes of a conventional interferometric radio telescope.
Like traditional radio dishes, these stations provide collecting area and raw
sensitivity as well as pointing and tracking capabilities. Fundamentally,
LOFAR differs from high-frequency radio telescopes in the fact that LOFAR
stations do not physically move. LOFAR operates by pointing and tracking using
combined signals from the individual antennas to form a phased array using a
combination of analog and digital beam forming techniques and in turn makes
LOFAR more flexible and agile. Station-level beam-forming allows for rapid re-
pointing of the telescope as well as the potential for multiple, simultaneous
observations from a given station. The resulting digitized, beam-formed data
from the stations can then be streamed to the central processing facility and
correlated to produce visibility’s for imaging applications and observation
analysis. In this research we use the processed observations generated by
I-LOFAR as seen in figure 1, to create a test set in which YOLOv2 can be
evaluated on.

Figure 1: I-LOFAR station located in Birr, Co.Offaly, Ireland. This an
international station consisting of 96 Low band antennas (right) operating
between 10-90MHz and 96 High band antennas (left) operating between
110-250MHz.

###  II-B Type III solar radio bursts

The Sun is an active star that produces large-scale events such as Coronal
Mass Ejections (CMEs) and solar flares. Radio emission is often associated
with these events in the form of radio bursts. These bursts are classified
into five main Types. Type I bursts are short duration narrowband bursts
associated with active regions. Type II bursts are slow frequency drifting
radio emissions thought to be excited by shock waves travelling through the
solar corona and they are associated with CMEs. Type III radio bursts are
rapid frequency drifting bursts which can sometimes be followed by continuum
emissions; these emissions are called Type V radio bursts. Type IV bursts are
broad continuum emissions with rapidly varying time structures. Solar radio
bursts are most often observed in dynamic spectra of frequency versus time.

Figure 2: Examples of real Type III solar radio bursts as seen in a dynamic
spectrum at a frequency range of 20-90MHz. This is used to in the test set to
evaluate YOLOv2s performance.

The most commonly occurring radio burst is the Type III, which is generally
short lasting a couple of seconds and structurally represent a vertical bright
strip in dynamic spectra at a frequency range of 10-100MHz, see Figure 2.
Although the structure of the Type III seems to be very basic, identifying
them is a complex task as they come in a variety of forms within the dynamic
spectra. For example, the Type III may be smooth or patchy, weak or strong,
superimposed on other radio bursts, standalone or in groups, or may be
embedded in strong radio frequency interference (RFI).

###  II-C Dataset

In order to train a useable model with YOLOv2, a large training dataset is
required. For our experiments, this dataset consisted of simulated Type III
data created using parametric models. Using these parametric models, we
produce Type III radio bursts that are random in number, grouping, intensity,
drift rate, and homogeneity. They also allow random variations in the
frequency range and duration’s of the bursts. We embed the bursts in a
background of simulated and random RFI channels, an example which can be found
in Figure 3.

Figure 3: A comparison between simulated Type III radio bursts (left) and
actual Type III radio bursts (right). The simulated Type III bursts are used
to train the deep learning model YOLOv2.

While producing these simulated Type III images we needed to label the
bounding box coordinates where the simulated Type III was randomly fixed into
the image. Using these simulated Type IIIs, a training set of 80,000 simulated
Type III SRBs was created to train the YOLOv2 model.

###  II-D Model Configuration

Once the dataset was created, we need to configure the model. We used a
framework called darkflow[17]; a Tensorflow implementation of darknet to
develop the model in YOLOv2.

TABLE I: The CNN architecture of YOLOv2. The CNN in YOLOv2 is altered in the
fully connected layers of the CNN. The fully connected layers are removed and
the detection and classifications are done by K-means classification for
improved accuracy. Darknet-19 Architecture
---
Type |  Filters |  Size/Stride |  Output
Convolutional |  32 |  3 X 3 |  224 X 224
Maxpool |  |  2 X 2/2 |  112 X 112
Convolutional |  64 |  3 X 3 |  112 X 112
Maxpool |  |  2 X 2/2 |  56 X 56
Convolutional |  128 |  3 X 3 |  56 X 56
Convolutional |  64 |  1 X 1 |  56 X 56
Convolutional |  128 |  3 X 3 |  56 X 56
Maxpool |  |  2 X 2/2 |  28 X 28
Convolutional |  256 |  3 X 3 |  28 X 28
Convolutional |  128 |  1 X 1 |  28 X 28
Convolutional |  256 |  3 X 3 |  28 X 28
Maxpool |  |  2 X 2/2 |  14 X 14
Convolutional |  512 |  3 X 3 |  14 X 14
Convolutional |  256 |  1 X 1 |  14 X 14
Convolutional |  512 |  3 X 3 |  14 X 14
Convolutional |  256 |  1 X 1 |  14 X 14
Convolutional |  512 |  3 X 3 |  14 X 14
Maxpool |  |  2 X 2/2 |  7 X 7
Convolutional |  1024 |  3 X 3 |  7 X 7
Convolutional |  512 |  1 X 1 |  7 X 7
Convolutional |  1024 |  3 X 3 |  7 X 7
Convolutional |  512 |  1 X 1 |  7 X 7
Convolutional |  1024 |  3 X 3 |  7 X 7
Convolutional |  1000 |  3 X 3 |  7 X 7
Avgpool |  |  Global |  1000
Softmax |  |  |

In YOLOv2 there are 19 convolutional layers and 5 maxpool layers as seen in
table 1. We needed to change the number of filters in the last convolutional
layer as we are only looking to detect one class being Type IIIs. To do this,
we use the function (1)

$filters=bounding*(classes+coords)$ (1)

where $filters$ is the number of filters in the last convolutional layer,
$bounding$ is the number of bounding boxes per grid cell, $classes$ is the
number of classes being detected, and $coords$ is the number of coordinates in
each bounding box. Here $bounding=5$, $classes=1$, and $coords=5$. The total
number of filters needed $filters$ at the final convolutional layer is 5 *
(1+5)=30.
We then altered the sizes of the bounding boxes using the anchor values. We
decided to set the bounding box width to be very narrow, rational being Type
IIIs are short lasting in terms of times often a couple of seconds. We then
set the height of the bounding box to be the height of 10-90MHz similarly seen
in the training set and dynamic spectrum. It’s apparent, Type III intensity’s
drift so low it cannot be seen by the naked eye as seen in Figure 4.

Figure 4: An Example case of a Type III intensity fading to the point of being
unidentifiable by the naked eye. Once the colours have been inverted the drift
of intensity can be seen.

In order to achieve real time frame rates with YOLOv2, we needed to decrease
the number of bounding boxes predicted on each test image from 845 to 405. To
do this, we needed to change the input size from 416x416 to 288x288. YOLOv2’s
convolutional layers downsample the image by a factor of 32 so by using an
input image of 416/32 we get an output feature map of 13 × 13 or sectioning an
image into a 13x13 square grid. This multiplied by the number of bounding
boxes per grid cell gives us the number of predictions per image 5x13x13=845.
In our case, we are feeding in a 288x288 input image. As a result we get a
feature map of 9x9 with 5 bounding boxes per square grid giving us 5x9x9=405
predictions per image. At 288×288 it runs at a maximum 77 fps with accuracy
comparable to Fast R-CNN. This input size complements the dataset input into
the YOLOv2 model as only one class in greyscale format is being identified.

##  III Results

###  III-A Training

The YOLOv2 model is trained to detect and classify Type III SRBs. The training
set consists of 80,000 simulated Type III images that are random in number,
grouping, intensity, drift rate, in-homogeneity, start-end frequency and
start-end time. The training set is also simultaneously labelled as it’s
created, allowing us to create a high volume training set along with a list of
text files containing bounding box coordinates. These text files are then
translated into XML to fit the darkflows training set requirements. Once
prerequisites have been complete the training set can be fed into the CNN for
training. Darkflow, by default, only uses 1 GPU or CPU for training so we had
to add the ability for Nvidia Scalable Link Interface (SLI) support to the
framework. The model is trained with a learning rate of 0.001, a momentum of
0.9 and Leaky Rectified Linear Unit (ReLU) as an activation function during
training. Learning parameters are also updated until convergence using
Stochastic Gradient Descent (SGD).

This research has been performed on a machine comprising of 2 x SLI inter-
connected GPU Nvidia Geforce RTX 2080 Ti, using Ubuntu 20.4.2 LTS on an AMD
Ryzen Threadripper 1950x with 32GB of RAM. For the training configuration we
are use 90% of GPU capacity for 200 epochs at a batch size of 32. With this
configuration, training took 4 days with loss decreasing with every iteration
as seen in Figure 5.

Figure 5: A plot showing loss decreasing with every iteration of training.
Loss is another method of evaluating how well YOLOv2 models the training
dataset.

###  III-B Test Set

The test set for the model is an 12 hour observation made by I-LOFAR on the
10th of September 2017. The raw data is processed and converted into an image
and plotted frequency in MHz from 0 to 250 over the duration of the
observation (time) as seen in figure 6. This image is then converted to
greyscale to match the characteristics of the training set so that colour
isn’t an influence on the models predictions. We focus on the 10-90 MHz range
in the observations as this shows the Type IIIs most prominent attributes. The
observation is then cut into 10 minute chunks to create a test set of 1331
images, containing about 15,000 Type III solar radio bursts. Once we have our
test set images, we needed to annotate our ground truth bounding box values.
This was done by using LabelImg, an annotating tool used to label objects
within an image. Once a Type III is labelled it’s corresponding bounding box
coordinates is stored in an XML file that will be used to compare the ground
truth coordinates to the models predicted coordinates.

Figure 6: An I-LOFAR observation made on the 10th of September 2017. This raw
data generated by I-LOFAR has been processed and plotted on a frequency range
of 0-250MHz.

###  III-C Model Evaluation

The YOLOv2 model’s performance was measured using the test set of 1331 images
from I-LOFAR described previously. The unit in which we represent our models
performance result is the f1-score. The f1-score highlights the balance
between precision and recall with precision meaning how good the model
predicted the location of an object and recall measuring how good the model
found and located all objects. We use the following functions to calculate the
f1-score of the model.

$Precision=\frac{TP}{TP+FP}$ (2) $Recall=\frac{TP}{TP+FN}$ (3)
$f1-score=2*\frac{Precision*Recall}{Precision+Recall}$ (4)

True positive (TP) and false positive (FP) are obtained from intersection over
union (IoU) from tested data. IoU compares the predicted bounding box to the
ground truth bounding box. A prediction is classified as TP if the IoU is
greater than 0.5, FP is if it’s less than 0.5. A False Negative (FN) is
specified for those images where the model failed to detect a known Type III
object. One important factor to take into account when evaluating the model
performance is confidence threshold. Confidence threshold measures how
confident the model is at predicting a certain object, in this case Type III
SRBs. The lower the confidence, the more detections made on a test image but
also the more false detections being made. After experimenting with different
thresholds it was found that the model was at it’s most optimized with
confidence threshold set to 0.1, see table 2. With this configuration, the
resulting f1-score is 82.63% for detecting Type III solar radio bursts.

TABLE II: From table 2 we can see that if we use higher threshold the
resulting accuracy will be lower. But if we lower the threshold to allow more
detections per test image the accuracy will be higher. However, in doing this
we can see an increase in false positives also. Model accuracy per confidence
threshold
---
Conf. thresh |  Precision |  Recall |  f1-score |  True Positive |  False Positive
0.5 |  10.79% |  98.35% |  19.44% |  1310 |  22
0.4 |  26.22% |  98.36% |  41.41% |  3185 |  53
0.3 |  46.33% |  98.13% |  63.04% |  5628 |  107
0.25 |  56.81% |  97.57% |  63.04% |  6910 |  172
0.2 |  68.20% |  96.75% |  80.01% |  8313 |  279
0.15 |  75.13% |  90.73% |  82.30% |  9220 |  942
0.1 |  83.42% |  81.85% |  82.63% |  9796 |  1947

## Conclusion

In this paper, we have shown that YOLOv2 is very good at detecting and
classifying Type III SRBs at high frame rates. This particular configuration
of YOLOv2 can achieve an accuracy of 82.63% on a real data observation
consisting of almost 15,000 Type III solar radio burst examples while also
achieving real-time frame rates (maximum 77 fps). This is a significant step
towards having an automated solar instrument that offers real time analysis.
While current state of the art algorithms such as Hough and randon transform
and CRAF offer excellent accuracy, these algorithms don’t have the benefit of
providing both accuracy and real time performance. Key to attaining the high
accuracy with the YOLO model is the quality of the training dataset which is
made up entirely of simulated Type III SRBs. We intend to create a new
training dataset with real Type III observations included with the
simulations. This will be used to test the robustness of the training set.
Issues associated with robustness would be the exclusion of embedded RFI and
general background noise as well as compensating for intensity fluctuations in
real data. In conclusion, we have shown that accurate real-time classification
of Type III SRBs is readily attainable.

## References

  * [1] R. P. Lin, “Energy release and particle acceleration in flares: Summary and future prospects,” Space Science Reviews, vol. 159, no. 1, p. 421, 2011\.
  * [2] V. N. Pick M, “Sixty five years of solar radio astronomy: flares, coronal mass ejections and sun-earth connection,” The Astronomy and Astrophysics Review, 2009.
  * [3] E. P. Carley, C. Baldovin, P. Benthem, M. M. Bisi, R. A. Fallows, P. T. Gallagher, M. Olberg, H. Rothkaehl, R. Vermeulen, N. Vilmer, and et al., “Radio observatories and instrumentation used in space weather science and operations,” Journal of Space Weather and Space Climate, vol. 10, 2020\.
  * [4] V. V. Lobzin, I. H. Cairns, and A. Zaslavsky, “Automatic recognition of type iii solar radio bursts in stereo/waves data for onboard real-time and archived data processing,” Journal of Geophysical Research A: Space Physics, vol. 119, no. 2, pp. 742––750, 2014.
  * [5] H. Salmane, R. Weber, K. Abed-Meraim, K. L. Klein, and X. Bonnin, “A method for the automated detection of solar radio bursts in dynamic spectra,” Journal of Space Weather and Space Climate, vol. 8, p. 0–18, 2018.
  * [6] L. Ma, Z. Chen, L. Xu, and Y. Yan, “Multimodal deep learning for solar radio burst classification,” Pattern Recognition, vol. 61, p. 573–582, 2017\.
  * [7] Y. C. Hou, Q. M. Zhang, S. W. Feng, Q. F. Du, C. L. Gao, Y. L. Zhao, and Q. Miao, “Identification and extraction of solar radio spikes based on deep learning,” Solar Physics, vol. 295, no. 10, 2020.
  * [8] Y. G. Zhang, V. Gajjar, G. Foster, A. Siemion, J. Cordes, C. Law, and Y. Wang, “Fast radio burst 121102 pulse detection and periodicity: A machine learning approach,” arXiv, 2018.
  * [9] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 2016-Decem, p. 770–778, 2016\.
  * [10] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 2016-December, p. 779–788, 2016.
  * [11] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C. Y. Fu, and A. C. Berg, “Ssd: Single shot multibox detector,” Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 9905 LNCS, p. 21–37, 2016.
  * [12] R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Region-based convolutional networks for accurate object detection and segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 38, no. 1, p. 142–158, 2016.
  * [13] R. Girshick, “Fast r-cnn,” Proceedings of the IEEE International Conference on Computer Vision, vol. 2015 Inter, p. 1440–1448, 2015.
  * [14] S. Ren, K. He, R. Girshick, and J. Sun, “Faster r-cnn: Towards real-time object detection with region proposal networks,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 39, no. 6, p. 1137–1149, 2017\.
  * [15] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask r-cnn,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 42, no. 2, p. 386–397, 2020.
  * [16] J. Redmon and A. Farhadi, “Yolo9000: Better, faster, stronger,” Proceedings - 30th IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, vol. 2017-Janua, p. 6517–6525, 2017.
  * [17] T. H. Trieu, “Darkflow,” GitHub Repository. Available online: https://github. com/thtrieu/darkflow (accessed on 14 February 2019), 2018.


11institutetext: Department of Information Science and Technology, Kyushu
University, Japan 11email<EMAIL_ADDRESS>
22institutetext: Department of Informatics, Kyushu University, Japan 22email:
<EMAIL_ADDRESS>
33institutetext: M&D Data Science Center, Tokyo Medical and Dental University,
Japan 33email<EMAIL_ADDRESS>
44institutetext: Department of Computer and Network Engineering, University of
Electro-Communications, Japan 44email<EMAIL_ADDRESS>
# Online algorithms for finding distinct substrings with length and multiple
prefix and suffix conditions

Laurentius Leonard 11 0000-0001-8477-7033
Shunsuke Inenaga 22 0000-0002-1833-010X
Hideo Bannai 33 0000-0002-6856-5185
Takuya Mieno 44 0000-0003-2922-9434

###### Abstract

Let two static sequences of strings $P$ and $S$, representing prefix and
suffix conditions respectively, be given as input for preprocessing. For the
query, let two positive integers $k_{1}$ and $k_{2}$ be given, as well as a
string $T$ given in an online manner, such that $T_{i}$ represents the
length-$i$ prefix of $T$ for $1\leq i\leq|T|$. In this paper we are interested
in computing the set $\mathit{ans_{i}}$ of distinct substrings $w$ of $T_{i}$
such that $k_{1}\leq|w|\leq k_{2}$, and $w$ contains some $p\in P$ as a prefix
and some $s\in S$ as a suffix. More specifically, the counting problem is to
output $|\mathit{ans_{i}}|$, whereas the reporting problem is to output all
elements of $\mathit{ans_{i}}$, for each iteration $i$. Let $\sigma$ denote
the alphabet size, and for a sequence of strings $A$, $\|A\|=\sum_{u\in
A}|u|$. Then, we show that after $O((\|P\|+\|S\|)\log\sigma)$-time
preprocessing, the solutions for the counting and reporting problems for each
iteration up to $i$ can be output in $O(|T_{i}|\log\sigma)$ and
$O(|T_{i}|\log\sigma+|\mathit{ans_{i}}|)$ total time. The preprocessing time
can be reduced to $O(\|P\|+\|S\|)$ for integer alphabets of size polynomial
with regard to $\|P\|+\|S\|$. Our algorithms have possible applications to
network traffic classification.

###### Keywords:

pattern matching counting algorithm suffix array suffix tree.

##  1 Introduction

Pattern matching has long been a central topic in the field of string
algorithms [4], leading to various applications, including DNA analysis in
bioinformatics [7, 13] as well as packet classification [3, 6, 22, 16] and
anti-spam email filtering in network security [17].

In this paper, we propose algorithms for counting and reporting distinct
substrings of an online text $T$ that have some $p\in P$ as a prefix and some
$s\in S$ as a suffix, and whose length is within the interval
$[k_{1}..k_{2}]$, where $P$ and $S$ are static sequences of strings given as
input for preprocessing, and integers $k_{1},k_{2}$ and the characters of $T$
are given as query. A similar yet different problem where patterns are given
in the form of a pair of a prefix and a suffix condition, i.e. $p\Sigma^{*}s$
patterns, rather than a pair of sequences where one is of prefixes and the
other of suffixes, is well-studied as the _followed-by_ problem [2, 14] or the
_Dictionary Recognition with One Gap (DROG)_ [1, 12, 21] problem. The problem
of this paper is also of importance with possible applications in network
traffic classification: All the application signatures in [20], for example,
can be expressed as an instance of our problem, as discussed in Section 4 and
demonstrated in Appendix 0.A. Note that while traffic classification via these
signatures was shown to be highly accurate, there are still cases of false
positives. When analyzing which patterns give false positives, we may be
interested in which patterns match the signatures, in which case the distinct
condition of our problem helps prevent wasting computation time on repeated
occurrences of each pattern. Another possible application is in computing the
distinct substrings of an online text $T$ whose lengths are at least $k$ and
belong to a $(k,r)$-TTSS language [19]; the words of length $\geq k$ in a
$(k,r)$-TTSS language defined by the 4-tuple $(I_{k},F_{k},T_{k,r},g)$ are the
words that have some element of $I_{k}$ as a prefix, some element of $F_{k}$
as a suffix, and includes, for each $t\in T_{k,r}$, at most $g(t)$ occurrences
of $t$ as substring. Here, the elements of $I_{k}$ and $F_{k}$ are strings of
length $k-1$ and the elements of $T_{k,r}$ are strings of length between 1 to
$k$ inclusive, and $g$ is a function that projects
$T_{k,r}\rightarrow\\{0,1,\cdots,r-1\\}$. A direct application of our
algorithms can check whether $w$ fulfills the prefix and suffix condition,
while the condition of restricted segments, i.e. the number of occurrences of
$T_{k,r}$ can also be considered by implementing the following modification:
for each iteration $i$, maintain the minimum start-index $\mathit{res}$ of the
suffix of $T_{i}$ that meets the condition of restricted segments, and use it
to exclude any suffixes longer than $T_{i}[\mathit{res}..i]$ from the
solution.

Our proposed algorithms take $O((\|P\|+\|S\|)\log\sigma)$ preprocessing time,
while processing $T$ itself in an online manner and outputting the solutions
up to iteration $i$ takes $O(|T_{i}|\log\sigma)$ and
$O(|T_{i}|\log\sigma+|\mathit{ans}_{i}|)$ cumulative time for the counting and
reporting problems respectively, using $O(|T_{i}|+\|P\|+\|S\|)$ working space.
Here, $T_{i}$ denotes the length-$i$ prefix of $T$, $\|P\|$ and $\|S\|$ denote
the total length of strings in $P$ and $S$ respectively, $\sigma$ is the
alphabet size, $|\mathit{ans}_{i}|$ is the number of substrings reported for
each $T_{i}$, and cumulative time refers to the total amount of running time
up to iteration $i$, as opposed to the running time of only iteration $i$. In
addition, the preprocessing time can be reduced to $O(\|P\|+\|S\|)$ in the
case of integer alphabets of size polynomial in $\|P\|+\|S\|$.

Also note that, the problems addressed in this paper differ from those of
[11], in which a different set of solution strings is output for each $p\in P$
where each solution must have that specific element of $P$ as a prefix, unlike
the problem in this paper where only a single set of solution strings is
output, where its elements can have any $p\in P$ as a prefix. Also, in [11]
there is only one suffix condition and no length condition, and the algorithm
is offline w.r.t. $T$, unlike in this paper.

##  2 Preliminaries and definitions

###  2.1 Strings

Let $\Sigma$ be an alphabet of size $\sigma$. An element of the set
$\Sigma^{*}$ is a string. The length of a string $w$ is denoted by $|w|$. The
empty string is denoted by $\varepsilon$. That is, $|\varepsilon|=0$. For a
string $w=pts$, $p$, $t$, and $s$ are called a _prefix_ , _substring_ , and
_suffix_ of $w$, respectively. A prefix $p$ (resp. suffix $s$) of a string $w$
is called a _proper prefix_ (resp. _proper suffix_) of $w$ if $|p|<|w|$ (resp.
$|s|<|w|$). For a string $w$, $w[i]$ denotes the $i$-th symbol of $w$ for
$1\leq i\leq|w|$, and $w[i..j]$ denotes the substring $w[i]w[i+1]\cdots w[j]$
for $1\leq i\leq j\leq|w|$. For a sequence $S$ of strings, let
$\|S\|=\sum_{u\in S}|u|$.

###  2.2 Suffix array and LCP array

The suffix array [15] of a string $w$ is a lexicographically sorted array of
suffixes of $w$, where each suffix is represented by its start-index. The LCP
array is an auxiliary array commonly used alongside the suffix array, that
stores the length of the longest common prefix of each adjacent pair of
suffixes in the suffix array. More specifically, if $\mathit{SA}$ and
$\mathit{LCP}$ are the suffix array and LCP array of the same string $w$, for
$x\in[2..|w|]$, $\mathit{LCP}[x]$ is the length of longest common prefix of
the suffixes $w[\mathit{SA}[x]..|w|]$ and $w[\mathit{SA}[x-1]..|w|]$. In this
paper, we will use suffix arrays for some strings, in which we denote by
$\mathit{LCP}$ the LCP array of the same string of the suffix array being
discussed. It is well-known that the suffix array and LCP array of a string
$w$ can be built in $O(|w|)$ time for integer alphabets of polynomial size in
$|w|$ [8, 9, 10], and in $O(|w|\log\sigma)$ time for general ordered alphabets
[23].

###  2.3 The problems

The problems considered in the paper are as follows.

######  Definition 1 (Online substring counting and reporting problem with
distinctness, multiple prefixes, multiple suffixes and length range
conditions)

Given two sequences of strings $P=(p_{1},\cdots,p_{n})$ and
$S=(s_{1},\cdots,s_{m})$, two integers $k_{1}$ and $k_{2}$, and a string $T$
given in an online manner (i.e., $T_{0}=\varepsilon$ and for each iteration
$i=1,\dots,|T|$, the $i$-th character is appended to $T_{i-1}$ to form
$T_{i}$), let $\mathit{ans}_{i}$ denote the set of distinct substrings of
$T_{i}$ that have some $p\in P$ as a prefix, some $s\in S$ as a suffix, and
whose length falls within the interval $[k_{1}..k_{2}]$.

  * The counting problem. On each iteration $i$, output $|\mathit{ans}_{i}|$.

  * The reporting problem. On each iteration $i$, output $\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$.

This paper excludes the empty string $\varepsilon$ from the solutions.

##  3 Algorithm

###  3.1 Sketch of algorithm

In this section, we describe the general idea of our algorithm. During each
iteration $i$, we need either to compute the size of
$\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$ to add it to the counting
solution, or to report all its elements. All elements of
$\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$ must be suffixes of $T_{i}$, and
thus for the suffix condition, $T_{i}$ itself must have some element of $S$ as
a suffix; otherwise clearly
$\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}=\emptyset$ and there is no need
to output a solution for the current iteration. Thus, let us consider the case
where $T_{i}$ has at least one element of $S$ as a suffix, and call the
shortest of them $s$. To help keep track of which suffixes of $T_{i}$ fulfill
the prefix condition, let us maintain a linked list $\mathit{pList}$ that
contains in increasing order, all distinct indices $j$ such that there is an
element of $P$ that occurs in $T_{i}$ with start-index $j$. Clearly, the
elements of $\mathit{pList}$ represent a bijection to the suffixes of $T_{i}$
that have an element of $P$ as a prefix, which are candidates for elements of
the solution set $\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$. Specifically,
for all $j\in\mathit{pList}$, the string
$u=T_{i}[j..i]\in\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$ iff $u$ fulfills
all the following conditions:

  1. (a)

$u$ has $s$ as a suffix.

  2. (b)

$u$ does not occur in $T_{i-1}$.

  3. (c)

$|u|\geq k_{1}$ .

  4. (d)

$|u|\leq k_{2}$.

Here, (a) $u$ has $s$ as a suffix iff $|u|\geq|s|$, and (b) $u$ does not occur
in $T_{i-1}$ iff $|u|>|\mathit{lrs}_{i}|$, where $\mathit{lrs}_{i}$ denotes
the longest repeating suffix of $T_{i}$, i.e. the longest suffix of $T_{i}$
that occurs at least twice in $T_{i}$. Thus, conditions (a) to (c) set a lower
bound for the length of suffixes of $T_{i}$ corresponding to elements of
$\mathit{pList}$ that can be a solution while condition (d) sets an upper
bound. If we visualize $\mathit{pList}$ horizontally as shown in Figure 1,
conditions (a) to (c) exclude some elements from the right while condition (d)
excludes some element from the left.

Figure 1: A visualized example of $\mathit{pList}$.

Take the maximum among the number of elements excluded by conditions (a) to
(c) and denote it by $\mathit{excludeRight}$, and denote the number of
elements excluded by condition (d) by $\mathit{excludeLeft}$. Then,
$|\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}|$ $=$
$\max{}(0,|\mathit{pList}|-\mathit{excludeLeft}-\mathit{excludeRight})$,
giving us the solution for the counting problem. For the reporting problem,
let us maintain $\mathit{start}$, a pointer to the smallest element in
$\mathit{pList}$ not excluded by condition (d). Then report all elements
traversed by starting at $\mathit{start}$ and moving to the right
$|\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}|-1$ times.

###### Example 1

The example in Figure 1 occurs when $T_{i}=\mathtt{coldcocoaold}$,
$P=(\mathtt{cave},\mathtt{coco},\mathtt{cocoa},\mathtt{d},\mathtt{oao},\mathtt{old})$,
$S=(\mathtt{aold},\mathtt{oaold})$, $k_{1}=3$, and $k_{2}=8$. Then,
$\mathit{pList}=(2,4,5,8,10,12)$. $10,12$ correspond to
$\mathtt{old},\mathtt{d}$ which are shorter than $s=\mathtt{aold}$, and thus
excluded by condition (a). $\mathit{lrs}_{i}=\mathtt{old}$ and thus condition
(b) also excludes $10$ and $12$, while condition (c) excludes only $12$ which
corresponds to $\mathtt{d}$. Therefore,
$\mathit{excludeRight}=\max{(2,2,1)}=2$. Meanwhile, condition (d) excludes $2$
and $4$ which correspond to $\mathtt{oldcocoaold}$ and $\mathtt{dcocoaold}$
and so $\mathit{excludeLeft}=2$.Thus, we have that
$|\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}|=\max{(0,6-2-2)}=2$. For the
reporting solution, we have that $\mathit{start}$ points to $5$. Traversing
$2$ elements starting from $8$ gives us $5$ and $8$, each corresponding to
$\mathtt{cocoaold}$ and $\mathtt{oaold}$, exactly the elements of
$\mathit{ans}_{i}\setminus\mathit{ans}_{i-1}$.

###  3.2 Removing redundant elements

We say that $p_{k}\in P$ is _redundant_ iff there exists $p_{k^{\prime}}$ of
$P$ s.t. either $p_{k}$ has $p_{k^{\prime}}$ as a proper prefix, or
$p_{k}=p_{k^{\prime}}\wedge k>k^{\prime}$. Similarly, $s_{k}\in S$ is
_redundant_ iff there exists $s_{k^{\prime}}$ of $S$ s.t. either $s_{k}$ has
$s_{k^{\prime}}$ as a proper suffix, or $s_{k}=s_{k^{\prime}}\wedge
k>k^{\prime}$.

It is not hard to see why they are called redundant; when multiple copies of
the same string exist in $P$, keeping only one copy suffices, and when
$p_{k}\in P$ has $p_{k^{\prime}}\in P$ as a proper prefix, the strings that
have $p_{k}$ as a prefix is a subset of strings that have $p_{k^{\prime}}$ as
a prefix, and thus the solution remains the same even if we delete $p_{k}$
from $P$. The same can be said for redundant elements of $S$.

As one part of the preprocessing, we rebuild $P$ and $S$ so that the redundant
elements are deleted. First, we describe how to rebuild $P$. Let
$\mathit{P_{seq}}=\$p_{1}\$p_{2}\$\cdots\$p_{n}\$$, where $\$\not\in\Sigma$
and $\$\prec c$ for all $c\in\Sigma$, and let $\mathit{SA}$ be the suffix
array of $\mathit{P_{seq}}$.

Then, each $\mathit{SA}[x]$ for $x\in[2..n+1]$ corresponds to the start-index
of $\$p$ in $\mathit{P_{seq}}$, for some $p\in P$. Starting from $x=2$, output
the corresponding $p$, namely the unique $p\in P$ s.t. $\$p\$$ occurs on index
$\mathit{SA}[x]$. Then, increment $x$ (at least once) until we have that
$\mathit{LCP}[x]<|\$p|$, i.e. until we find an index that corresponds to
$\$p^{\prime}$ where $p^{\prime}\in P$ does not have $p$ as a prefix. Output
the element of $P$ corresponding to the new $x$, then again increment $x$ in
the same manner. Repeat this until $x>n+1$, and we have that all non-redundant
elements of $P$ are output.

Other than the construction of the suffix array and LCP array, clearly this
takes $O(\|P\|)$ time, and the same method can be used to compute non-
redundant elements of $S$: Let $S^{-1}=(s_{1}^{-1},\ldots,s_{m}^{-1})$ be the
sequence of reversed elements of $S$, then apply the above algorithm to
$S^{-1}$ and reverse each string in the output to get the non-redundant
elements of $S$. Thus, both $P$ and $S$ are rebuilt to exclude redundant
elements in $O(\|P\|+\|S\|)$ time, in addition to the construction time of the
suffix array and LCP array, which depends on the alphabet.

###### Example 2

Let $P=(\mathtt{abc},\mathtt{ab},\mathtt{acc},\mathtt{ab},\mathtt{cab})$.
Then, $\mathit{P_{seq}}=\mathtt{\$abc\$ab\$acc\$ab\$cab\$}$ and we have the
table as shown in Figure 2. $x=2$ corresponds to the occurrence of
$\mathtt{\$ab}$, which occurs on $\mathit{P_{seq}}$ on index
$\mathit{SA}[2]=5$. Thus, $\mathtt{ab}$ is determined to be non-redundant, and
we have that $|\mathtt{\$ab}|=3$, so increment $x$ until we have that
$\mathit{LCP}[x]<3$. This skips over $x=3,4$, correctly determining their
corresponding elements of $P$, namely the second $\mathtt{ab}$ and
$\mathtt{abc}$ to be redundant. When $x=5$, $\mathit{LCP}[x]=2<3$ and so the
corresponding $p=\mathtt{acc}$ is output. Similarly, for $x=6$,
$\mathit{LCP}[x]=1<4$ and thus $\mathtt{cab}$ is output. Afterwards, $x$ is
incremented beyond the interval $[2..n+1]$ and thus the algorithm terminates
and the non-redundant elements $(\mathtt{ab},\mathtt{acc},\mathtt{cab})$ are
output.

Figure 2: The table for $\mathit{P_{seq}}=\$abc\$ab\$acc\$ab\$cab$

For the rest of the paper, we will assume that the above preprocessing is done
and thus $P=(p_{1},\ldots,p_{n})$ and $S=(s_{1},\ldots,s_{m})$ from this point
refer to the rebuilt sequences that have no redundant elements.

###  3.3 Detecting $P$ and $S$ occurrences

As discussed, on each iteration $i$ we need to detect whenever an element of
$S$ occurs as a suffix of $T_{i}$. Additionally, we will also need to detect
when an element of $P$ occurs as a suffix of $T_{i}$, in order to maintain
$\mathit{pList}$. This can be done by building an Aho-Corasick automaton for
$P$ and $S$ separately. Constructing both automata takes
$O((\|P\|+\|S\|)\log\sigma)$ preprocessing time in general, and
$O(\|P\|+\|S\|)$ time in the case of integer alphabets of size polynomial with
regard to $\|P\|+\|S\|$ [5]. Running each automaton up to iteration $i$ takes
$O(|T_{i}|\log\sigma+\mathit{occ})$ cumulative time, where $\mathit{occ}$ is
the number of occurrences detected. Here, the occurrences of elements of $P$
in $T_{i}$ must have distinct start-indices, as two occurrences with a shared
start-index imply that one of them is redundant. Similarly, occurrences of
elements of $S$ must have distinct end-indices and thus $\mathit{occ}\in
O(|T_{i}|)$ for both automata, and so the cumulative running time becomes
$O(|T_{i}|\log\sigma)$.

###  3.4 Maintaining $\mathit{pList}$

To maintain $\mathit{pList}$, whenever some $p\in P$ is detected to occur as a
suffix of $T_{i}$, its start-index $j$ needs to be added to $\mathit{pList}$
while maintaining the increasing order. Doing this naively would take
$O(|\mathit{pList}|)=O(|T_{i}|)$ time for every insertion which gives
quadratic time overall, so a more efficient scheme is necessary.

During any iteration $i$, the elements of $P$ that occur as suffixes of
$T_{i}$ are detected in decreasing order of length, because we use Aho-
Corasick automaton. Let $p$ be such an element detected, and $j$ be the start-
index of its occurrence, i.e. $j=i-|p|+1$. We need to add $j$ into
$\mathit{pList}$ so that the increasing order of its elements are maintained.

To do that, we need to find the minimum $j^{\prime}$ among the current
elements of $\mathit{pList}$ such that $j^{\prime}>j$. Here, $j^{\prime}$
being an element of $\mathit{pList}$ implies that $j^{\prime}$ corresponds to
an occurrence of $p^{\prime}\in P$ starting at $j^{\prime}$ and ending at some
$i^{\prime}\leq i$. We can see that in fact $i^{\prime}<i$, for if
$i^{\prime}=i$, $j^{\prime}$ was added to $\mathit{pList}$ in the current
iteration $i$ before $j$, while $j^{\prime}>j\wedge i^{\prime}=i$ implies
$|p^{\prime}|<|p|$, contradicting the fact that the Aho-Corasick automaton
detects the occurrences in decreasing order of length. Thus, $j^{\prime}>j$
and $i^{\prime}<i$, meaning the occurrence of $p^{\prime}$ falls completely
within $p[2..|p|-1]$.

Our scheme is then as follows: We precompute, for each $p_{k}\in P$, the
minimum value $y$ such that there is some $p_{k^{\prime}}\in P$ that occurs in
$p_{k}[2..|p_{k}|-1]$ on start-index $y+1$. If there is no such
$p_{k^{\prime}}$, then let $y=\infty$. Then, we will store the values on the
array $\mathit{successorOffset}$ that maps each $p_{k}\in P$ to its
corresponding $y$.

Additionally, maintain also an array $\mathit{pArray}$ such that
$\mathit{pArray}[j]$ points to the element of $\mathit{pList}$ whose value is
$j$ if it exists, or $\mathit{null}$ otherwise. Then, whenever some $p_{k}\in
P$ occurs as a suffix of $T_{i}$ with start-index $j$, we can just add $j$
into $\mathit{pList}$ exactly before the element pointed to by
$pArray[j+\mathit{successorOffset}[k]]$. Clearly, once
$\mathit{successorOffset}$ is computed, adding each element of
$\mathit{pList}$ takes only constant time and thus maintaining
$\mathit{pList}$ and $\mathit{pArray}$ takes cumulative $O(|T_{i}|)$ time, as
the number of elements of $\mathit{pList}$ for any given iteration $i$ is
bounded by $|T_{i}|$.

##### Computing $\mathit{successorOffset}$.

Let $\mathit{P_{concat}}=p_{1}\$p_{2}\$\cdots\$p_{n}\$$. Note that it differs
from $\mathit{P_{seq}}$ not only with regard to the positioning of $\$$, but
also in that redundant elements of $P$ are not included. For $k\in[1..n]$, let
$\mathit{PL}[k]$ denote the start-index of $p_{k}$ in $\mathit{P_{concat}}$,
i.e. $\mathit{PL}[k]=1+\sum_{k^{\prime}\in[1..k)}{(|p_{k^{\prime}}|+1)}$. For
$j\in[1..|\mathit{P_{concat}}|]$, let $\mathit{PI}[j]$ be the index of the
element of $P$ covering index $j$ in $\mathit{P_{concat}}$. That is,
$\mathit{PI}[j]=k$ where $j\in[\mathit{PL}[k]..\mathit{PL}[k]+|p_{k}|-1]$ if
such $k$ exists, otherwise $\mathit{PI}[j]=\mathit{null}$. See Figure 3 for an
example.

Figure 3: Example of $\mathit{PL}$ and $\mathit{PI}$.

Next, we construct the suffix array $\mathit{SAP}$ of $\mathit{P_{concat}}$.
We can then compute $\mathit{successorOffset}$. The general idea is that for
each $p_{k}$, we find all its occurrences using the suffix array, and when the
occurrence falls within some $p_{k^{\prime}}\in P$, we update
$\mathit{successorOffset}[k^{\prime}]$. A more detailed description is as
follows:

  * •

Initialize $\mathit{successorOffset}[k]=\infty$ for all $p_{k}\in P$ .

  * •

Using the LCP array, find the subinterval $[\ell..r]$ in $\mathit{SAP}$ that
corresponds to occurrences of $p_{k}$, i.e. each of $\mathit{SAP}[x]$ for all
$x\in[\ell..r]$ corresponds to some start-index of occurrences of $p_{k}$ in
$\mathit{P_{concat}}$.

  * •

For each such occurrence, whenever it falls within
$p_{k^{\prime}}[2..|p_{k^{\prime}}|-1]$ for some $p_{k^{\prime}}\in P$, then
assign to $\mathit{successorOffset}[k^{\prime}]$ the minimum value between
itself and the offset distance $\mathit{SAP}[x]-\mathit{PL}[k^{\prime}]$.
Formally, for all $x\in[\ell..r]$, if
$\mathit{P_{concat}}[\mathit{SAP}[x]+|p|]\neq\$$, assign to
$\mathit{successorOffset}[\mathit{PI}[\mathit{SAP}[x]]]$ the following value:

$\min{(\mathit{successorOffset}[\mathit{PI}[\mathit{SAP}[x]]],\mathit{SAP}[x]-\mathit{PL}[\mathit{PI}[\mathit{SAP}[x]]])}$
(1)

Computing $\mathit{P_{concat}}$ and $\mathit{PI}$ trivially takes $O(\|P\|)$
time. To compute the subinterval corresponding to occurrences of $p_{k}$,
simply find the longest subinterval $[\ell..r]$ of $\mathit{SAP}$ that
includes the index $\mathit{SAP}^{-1}[\mathit{PL}[k]]$ and
$\mathit{LCP}[x]\geq|p_{k}|$ for all $x\in[\ell+1..r]$. This takes linear time
w.r.t. to the subinterval length, which all adds up to the total number of
occurrences of elements of $P$ in $\mathit{P_{concat}}$. As no occurrence may
share a start-index, it is bounded by $|\mathit{P_{concat}}|\in O(\|P\|)$.
Thus, this preprocessing takes $O(\|P\|)$ time, in addition to the time
required to compute $\mathit{SAP}$ and $\mathit{LCP}$ which depends on the
alphabet.

###  3.5 Computing $\mathit{excludeLeft}$, $\mathit{excludeRight}$, and
$\mathit{start}$

As discussed, the algorithm runs the Aho-Corasick automaton for $S$ to check
whether there is some $s\in S$ that occurs as a suffix of $T_{i}$, for each
iteration $i$. In case such $s$ exists, all of the computations below are
performed, otherwise only maintaining $\mathit{start}$ is necessary.

####  3.5.1 $\mathit{excludeRight}$

##### Exclusion by $s$.

When $s\in S$ occurs as a suffix of $T_{i}$, the suffix $u=T_{i}[j..i]$ for
each $j\in\mathit{pList}$ is excluded from the solution iff $j>i-|s|+1$. That
is, $u$ fails to meet condition (a) of Section 3.1 iff $j$ corresponds to an
occurrence of $p\in P$ in $s[2..|s|]$. Thus, if we preprocess the number of
occurrences of elements of $P$ that occur within $s[2..|s|]$ for each $s\in
S$, we can compute the number of elements of $\mathit{pList}$ excluded by $s$
in constant time.

###### Example 3

In the case shown by Example 1, the elements 10 and 12 are excluded by $s$,
which we can obtain from the fact that there are two occurrences of elements
of $P$ within $s[2..|s|]=\mathtt{old}$; one each of $\mathtt{old}$ and
$\mathtt{d}$.

A suffix array-based approach that is similar to what we used to compute
$\mathit{successorOffset}$ can be used here. Let

$\mathit{PS_{concat}}=p_{1}\$\cdots p_{n}\$s_{1}[2..|s_{1}|]\$\cdots
s_{m}[2..|s_{m}|],$

and construct its suffix array $\mathit{SAPS}$. Define an array $\mathit{SI}$
such that $\mathit{SI}[j]=k$ when $j$ belongs to the part made up by
$s_{k}[2..|s_{k}|]$ in $\mathit{PS_{concat}}$, similar to $\mathit{PI}$. Then,
for all $p\in P$, compute the subinterval $[\ell..r]$ in $\mathit{SAPS}$
corresponding to suffixes of $\mathit{PS_{concat}}$ that start with $p$, again
using $\mathit{SAPS}^{-1}$ and LCP arrays. For all $x\in[\ell..r]$, then
increment $\mathit{sPCount}[\mathit{SI}[\mathit{SAPS}[x]]]$ by one, where
$\mathit{sPCount}$ is an array that maps each $s\in S$ to the number of
occurrences of elements of $P$ that occur within $s[2..|s|]$. Naturally,
$\mathit{sPCount}$ initially maps all elements of $S$ to zero before the
counts are incremented.

The total of size of subintervals is bounded by the number of occurrences of
elements of $P$ in $\mathit{PS_{concat}}$, which is
$O(|\mathit{PS_{concat}}|)=O(\|P\|+\|S\|)$. Thus, computing
$\mathit{successorOffset}$ takes $O(\|P\|+\|S\|)$ time, in addition to the
construction time of $\mathit{SAPS}$ and $\mathit{LCP}$ which depends on the
alphabet. After the preprocessing, the number of elements of $\mathit{pList}$
excluded by $s$ when $s$ occurs as a suffix of $T_{i}$ can be computed in
constant time by simply referring to $\mathit{sPCount}[s]$.

##### Exclusion by $\mathit{lrs}_{i}$.

It is known that Ukkonen’s algorithm[23] maintains the locus of
$\mathit{lrs}_{i}$ during each iteration $i$ where it builds the suffix tree
of $T_{i}$, and thus we can compute its start-index $i-|\mathit{lrs}_{i}|+1$
by simply running Ukkonen’s algorithm. Then, clearly $j\in\mathit{pList}$ is
excluded by condition (b) of Section 3.1 iff $j\geq i-|\mathit{lrs}_{i}|+1$.
Furthermore, the start-index is non-decreasing between iterations, i.e.
$i-|\mathit{lrs}_{i}|+1\geq(i-1)-|\mathit{lrs}_{i-1}|+1$ for any iteration
$i$. Thus, we can always maintain the count $|\\{j\in\mathit{pList}\mid
j\in[i-|\mathit{lrs}_{i}|+1..i]\\}|$ for each $i$ as follows:

  * •

During some iterations, $i-|\mathit{lrs}_{i}|+1>(i-1)-|\mathit{lrs}_{i-1}|+1$
which we will know from Ukkonen’s algorithm. In that case, for each
$j\in[(i-1)-|\mathit{lrs}_{i-1}|+2..i-|\mathit{lrs}_{i}|+1]$ such that
$j\in\mathit{pList}$, decrement the count by one.

  * •

Whenever a new element $j$ is added to $\mathit{pList}$ such that $j\in
i-|\mathit{lrs}_{i}|+1$, increment the count.

We can check whether $j\in\mathit{pList}$ in constant time for any $j$ using
$\mathit{pArray}$, and both the left end $i-|\mathit{lrs}_{i}|+1$ and right
end $i$ of the interval only ever increases and is within $1$ to $i$, so the
above method takes cumulative $O(|T_{i}|)$ time, dominated by the runtime of
Ukkonen’s algorithm which is cumulative $O(|T_{i}|\log\sigma)$ time.

##### Exclusion by $k_{1}$.

The same approach can be used to find the number of elements of
$\mathit{pList}$ excluded by $k_{1}$: we maintain the number of
$j\in\mathit{pList}$ such that $i-j+1<k_{1}\Leftrightarrow j\in[i-k_{1}+2..i]$
for each iteration $i$. Since both endpoints of this interval can only
increase and are always between $1$ to $i$ inclusive, maintaining the count
can be done in cumulative $O(|T_{i}|)$ time.

####  3.5.2 $\mathit{excludeLeft}$ and $\mathit{start}$

##### Maintaining $\mathit{excludeLeft}$.

Similarly, $\mathit{excludeLeft}$ is the number of $j\in\mathit{pList}$ such
that $i-j+1>k_{2}\Leftrightarrow j\in[1..i-k_{2}]$, and the same approach can
be used to maintain this count in $O(|T_{i}|)$ time.

##### Maintaining $\mathit{start}$.

We can easily maintain $\mathit{start}$ so that it points to the minimum
element of $\mathit{pList}$ of value at least $i+1-k_{2}$ in cumulative
$O(|T_{i}|)$ time as follows:

  * •

Initialize $\mathit{start}$ to $\mathit{null}$.

  * •

During each iteration $i$, if $\mathit{start}$ is not $\mathit{null}$ and
$\mathit{start}$ < $i+1-k_{2}$, then let $\mathit{start}$ point to the next
element in $\mathit{pList}$.

  * •

Whenever a new element $j$ is added to $\mathit{pList}$ such that $j\geq
i+1-k_{2}$, if $\mathit{start}=\mathit{null}$ or $\mathit{start}>j$, let
$\mathit{start}$ point to $j$.

###  3.6 Summarizing the algorithm

In the preprocessing, suffix arrays and LCP arrays are constructed and used to
remove redundant elements of $P$ and $S$, as well as compute
$\mathit{successorOffset}$. Additionally, Aho-Corasick automata for $P$ and
$S$ are built. For general ordered alphabets, this preprocessing takes
$O((\|P\|+\|S\|)\log\sigma)$ time, with the construction time for Aho-Corasick
automata and suffix arrays being the bottleneck. In the case of integer
alphabets of size polynomial w.r.t. $\|P\|+\|S\|$, the construction times, and
consequently the whole preprocessing time, can be reduced to $O(\|P\|+\|S\|)$
time.

For the query processing time up to any iteration $i$, $\mathit{excludeLeft}$,
$\mathit{excludeRight}$, and $\mathit{start}$ are computed in
$O(|T_{i}|\log\sigma)$ cumulative time, giving us the solution for the
counting problem. For the reporting problem, we traverse a total of
$|\mathit{ans}_{i}|$ elements in $\mathit{pList}$, giving us
$O(|T_{i}|\log\sigma+|\mathit{ans}_{i}|)$ cumulative time, assuming the
solution strings are output in the form of index pairs.

Additionally, all the data structures used require only
$O(|T_{i}|+\|P\|+\|S\|)$ total working space. Thus, we have the following
results.

###### Theorem 3.1

There are algorithms that solve the counting and solving problems from
Definition 1 for general ordered alphabets, such that after
$O((\|P\|+\|S\|)\log\sigma)$ preprocessing time, the solutions are output for
each iteration up to $i$ in $O(|T_{i}|\log\sigma)$ time for the counting
problem, and $O(|T_{i}|\log\sigma+|\mathit{ans}_{i}|)$ time for the reporting
problem, using $O(|T_{i}|+\|P\|+\|S\|)$ total working space.

###### Corollary 1

The preprocessing time in Theorem 3.1 can be reduced to $O(\|P\|+\|S\|)$ time
in case of integer alphabets of size polynomial with regard to $\|P\|+\|S\|$.

##  4 Applying the algorithm for traffic classification

We show in Appendix 0.A the input sets that match each of the application
signatures described in [20].

Note that the signatures shown in [20] are generally characterized in the
format of $p$ followed by $s$, where $p\in P$ and $s\in S$ for sets or lists
$P$ and $S$. In general, this differs from our problem in that occurrences of
$p$ and $s$ overlapping should not be counted as a match. For example,
$\mathtt{ab}$ followed by $\mathtt{bc}$ means $\mathtt{abc}$ should not be
counted as a match, while our algorithms do count this as a match.
Nevertheless, such matches, which would be erroneous in the context of
implementing the signatures, do not occur with the input sets listed in
Appendix 0.A, as the elements of $P$ and $S$ simply cannot overlap. For
example, with Gnutella signatures each element of $P$ ends with the character
$\mathtt{:}$, which no element of $S$ contains, so no string $w\in\Sigma^{*}$
exists such that has some $p\in P$ as prefix, $s\in S$ as suffix, and $p$ and
$s$ overlap (i.e. $|w|<|p|+|s|$). Note also that this problem also differs
from the followed-by problem of [2, 14], which does share the intolerance of
such overlaps, but differs in that the inputs are given as pairs of $p$ and
$s$ rather than pair of sets or lists $P$ and $S$. Naively solving the pair-
of-sets problem using algorithms for the pairs of $p$ and $s$ problem would
take $|P|\times|S|$ queries, as we need one query for each pair of $p\in P$
and $s\in S$. This is clearly inefficient for large $|P|$ and $|S|$, and hence
the necessity remains for our proposed algorithms.

##  5 Conclusion and future work

In this paper, we proposed online algorithms for counting and reporting all
distinct substrings of an online text $T$ that has some $p\in P$ as a prefix,
some $s\in S$ as a suffix, and whose length is within the interval
$[k_{1}..k_{2}]$, where $P$ and $S$ are static sequences of strings given as
input for preprocessing, while positive integers $k_{1},k_{2}$ and the
characters of $T$ are given as query. Our algorithms take
$O((\|P\|+\|S\|)\log\sigma)$ preprocessing time for general ordered alphabets,
which is reduced to $O(\|P\|+\|S\|)$ time for integer alphabets of size
polynomial w.r.t. $\|P\|+\|S\|$. The computation up to the $i$-th character of
$T$ takes $O(|T_{i}|\log\sigma)$ cumulative time for the counting problem, and
$O(|T_{i}|\log\sigma+|\mathit{ans}_{i}|)$ cumulative time for the reporting
problem. Furthermore, we have shown that it has possible applications in
traffic classification, by showing that all of the application signatures in
[20] can be represented as input sets of our proposed problems.

A few problems remain to be considered as future work:

  * •

As the discussion in Section 4 implies, solving the problem where the prefix
and suffix strings are not allowed to overlap, i.e. substrings are in the form
of $p\Sigma^{k}s$, where $k\in[k_{1}..k_{2}],p\in P,s\in S$, while retaining
the distinctness condition as well as that the input sets be given as pairs of
lists $P$ and $S$, can be useful in case we have a signature implemented with
$P,S$ such that there does exist a string $w$ such that $w$ has $p\in P$ as
prefix, $s\in S$ as suffix and $|w|<|p|+|s|$, and we want to exclude such $w$
from matches. Is it possible to devise an algorithm that solve this problem
efficiently?

  * •

In practice, how do the running times of our algorithms compare to the
signature implementations used in [20]?

## Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP20H04141 (HB) and
JP22H03551 (SI), and by JST PRESTO Grant Number JPMJPR1922 (SI).

## References

  * [1] Amir, A., Levy, A., Porat, E., Shalom, B.R.: Online recognition of dictionary with one gap. Information and Computation 275, 104633 (2020)
  * [2] Baeza-Yates, R.A., Gonnet, G.H.: Fast text searching for regular expressions or automaton searching on tries. Journal of the ACM (JACM) 43(6), 915–936 (1996)
  * [3] Choi, Y.H., Jung, M.Y., Seo, S.W.: L+ 1-mwm: A fast pattern matching algorithm for high-speed packet filtering. In: IEEE INFOCOM 2008-The 27th Conference on Computer Communications. pp. 2288–2296. IEEE (2008)
  * [4] Crochemore, M., Rytter, W.: Text algorithms. Maxime Crochemore (1994)
  * [5] Dori, S., Landau, G.M.: Construction of Aho Corasick automaton in linear time for integer alphabets. In: Apostolico, A., Crochemore, M., Park, K. (eds.) Combinatorial Pattern Matching. pp. 168–177. Springer Berlin Heidelberg, Berlin, Heidelberg (2005)
  * [6] Fuchino, T., Harada, T., Tanaka, K., Mikawa, K.: Acceleration of packet classification using adjacency list of rules. 2019 28th International Conference on Computer Communication and Networks (ICCCN) (7 2019). https://doi.org/10.1109/icccn.2019.8846923
  * [7] Gusfield, D.: Algorithms on Strings, Trees, and Sequences - Computer Science and Computational Biology. Cambridge University Press (1997). https://doi.org/10.1017/cbo9780511574931
  * [8] Kärkkäinen, J., Sanders, P., Burkhardt, S.: Linear work suffix array construction. Journal of the ACM (JACM) 53(6), 918–936 (2006)
  * [9] Kasai, T., Lee, G., Arimura, H., Arikawa, S., Park, K.: Linear-time longest-common-prefix computation in suffix arrays and its applications. In: Annual Symposium on Combinatorial Pattern Matching. pp. 181–192. Springer (2001)
  * [10] Kim, D.K., Sim, J.S., Park, H., Park, K.: Constructing suffix arrays in linear time. Journal of Discrete Algorithms 3(2), 126–142 (2005). https://doi.org/10.1016/j.jda.2004.08.019
  * [11] Leonard, L., Tanaka, K.: Suffix tree-based linear algorithms for multiple prefixes, single suffix counting and listing problems (2022). https://doi.org/10.48550/ARXIV.2203.16908
  * [12] Levy, A., Shalom, B.R.: Online parameterized dictionary matching with one gap. Theoretical Computer Science 845, 208–229 (2020). https://doi.org/10.1016/j.tcs.2020.09.016
  * [13] Makinen, V., Belazzougui, D., Cunial, F., Tomescu, A.I.: Genome-scale algorithm design. Cambridge University Press, Cambridge, England (May 2015)
  * [14] Manber, U., Baeza-Yates, R.: An algorithm for string matching with a sequence of don’t cares. Information Processing Letters 37(3), 133–136 (1991). https://doi.org/10.1016/0020-0190(91)90032-D
  * [15] Manber, U., Myers, G.: Suffix arrays: A new method for on-line string searches. SIAM Journal on Computing 22(5), 935–948 (1993). https://doi.org/10.1137/0222058
  * [16] Mikawa, K., Tanaka, K.: Run-based trie involving the structure of arbitrary bitmask rules. IEICE Transactions on Information and Systems E98.D(6), 1206–1212 (2015). https://doi.org/10.1587/transinf.2013EDP7087
  * [17] Pampapathi, R., Mirkin, B., Levene, M.: A suffix tree approach to anti-spam email filtering. Machine Learning 65(1), 309–338 (2006)
  * [18] Roesch, M., et al.: Snort: Lightweight intrusion detection for networks. In: Lisa. vol. 99, pp. 229–238 (1999)
  * [19] Ruiz, J., España, S., García, P.: Locally threshold testable languages in strict sense: Application to the inference problem. In: International Colloquium on Grammatical Inference. pp. 150–161. Springer (1998)
  * [20] Sen, S., Spatscheck, O., Wang, D.: Accurate, scalable in-network identification of p2p traffic using application signatures. In: Proceedings of the 13th International Conference on World Wide Web. p. 512–521. WWW ’04, Association for Computing Machinery, New York, NY, USA (2004). https://doi.org/10.1145/988672.988742
  * [21] Shalom, B.R.: Parameterized dictionary matching and recognition with one gap. Theoretical Computer Science 854, 1–16 (2021). https://doi.org/10.1016/j.tcs.2020.11.017
  * [22] Tongaonkar, A.S.: Fast pattern-matching techniques for packet filtering. Ph.D. thesis, Stony Brook University (2004)
  * [23] Ukkonen, E.: On-line construction of suffix trees. Algorithmica 14(3), 249–260 (Sep 1995). https://doi.org/10.1007/BF01206331

##  Appendix 0.A Appendix

Below, we show the input sets that match the each of the application
signatures described in [20].

Table 1: Input sets corresponding to application signatures Application |  List of $P$ elements |  List of $S$ elements |  $k_{1}$ |  $k_{2}$
---|---|---|---|---
Gnutella |  $\mathtt{User-Agent:},$ $\mathtt{UserAgent:},$ $\mathtt{Server:}$ |  $\mathtt{LimeWire},$ $\mathtt{BearShare},$ $\mathtt{Gnucleus},$ $\mathtt{MorpheusOS},$ $\mathtt{XoloX},$ $\mathtt{MorpheusPE},$ $\mathtt{gtkgnutella},$ $\mathtt{Acquisition},$ $\mathtt{Mutella-0.4.1},$ $\mathtt{MyNapster},$ $\mathtt{Mutella0.4.1},$ $\mathtt{Mutella-0.4},$ $\mathtt{Qtella},$ $\mathtt{AquaLime},$ $\mathtt{NapShare},$ $\mathtt{Comeback},$ $\mathtt{Go},$ $\mathtt{PHEX},$ $\mathtt{SwapNut},$ $\mathtt{Mutella-0.4.0},$ $\mathtt{Shareaza},$ $\mathtt{Mutella-0.3.9b},$ $\mathtt{Morpheus},$ $\mathtt{FreeWire},$ $\mathtt{Openext},$ $\mathtt{Mutella-0.3.3},$ $\mathtt{Phex}$ |  $1$ |  $\infty$
eDonkey |  $\mathtt{0xe3}$ (in hex) |  (the packet length) |  5-byte long
DirectConnect |  $\mathtt{\$MyNick},$ $\mathtt{\$Lock},$ $\mathtt{\$Key},$ $\mathtt{\$Direction},$ $\mathtt{\$GetListLen},$ $\mathtt{\$ListLen},$ $\mathtt{\$MaxedOut},$ $\mathtt{\$Error},$ $\mathtt{\$Send},$ $\mathtt{\$Get},$ $\mathtt{\$FileLength},$ $\mathtt{\$Canceled},$ $\mathtt{\$HubName},$ $\mathtt{\$ValidateNick},$ $\mathtt{\$ValidateDenide},$ $\mathtt{\$GetPass},$ $\mathtt{\$MyPass},$ $\mathtt{\$BadPass},$ $\mathtt{\$Version},$ $\mathtt{\$Hello},$ $\mathtt{\$LogedIn},$ $\mathtt{\$MyINFO},$ $\mathtt{\$GetINFO},$ $\mathtt{\$GetNickList},$ $\mathtt{\$NickList},$ $\mathtt{\$OpList},$ $\mathtt{\$To},$ $\mathtt{\$ConnectToMe},$ $\mathtt{\$MultiConnectToMe},$ $\mathtt{\$RevConnectToMe},$ $\mathtt{\$Search},$ $\mathtt{\$MultiSearch},$ $\mathtt{\$SR},$ $\mathtt{\$Kick},$ $\mathtt{\$OpForceMove},$ $\mathtt{\$ForceMove},$ $\mathtt{\$Quit}$ |  $\mathtt{|}$ |  $1$ |  $\infty$
BitTorrent |  the 20-byte string where the first byte is 19 (0x13) and the next 19 bytes are the string $\mathtt{19BitTorrent\;\,protocol}$ |  20-byte long
Kazaa |  $\mathtt{GET}$, $\mathtt{HTTP}$ |  $\mathtt{X-Kazaa}$ |  $1$ |  $\infty$