myequation

$\BODY$

(1)

#
Real-time Automated Answer Scoring

Akash Nagaraj Department of Computer Science PES University Bangalore, India
<EMAIL_ADDRESS>Mukund Sood Department of Computer Science PES University
Bangalore, India<EMAIL_ADDRESS>Dr. Gowri Srinivasa Department of
Computer Science PES University Bangalore, India<EMAIL_ADDRESS>
###### Abstract

In recent years, the role of big data analytics has exponentially grown, and
is now slowly making its way into the education industry. Several attempts are
being made in this sphere in order to improve the quality of education being
provided to students while simultaneously improving teaching techniques to
reduce the burden put on the teacher. While the implications of such a
collaboration have in fact been carried out before, automated scoring of
answers has been explored to a rather limited extent.

The biggest obstacle to choosing constructed-response assessments over the
traditional multiple-choice assessments is the large cost, and effort that
comes with their scoring, and this is precisely the issue that this project
aims to solve.The aim is to accept raw-input from the student in the form of
their answer, pre-process these answers if required, and automatically score
the answers. Evaluation of these answers will be based purely on the previous
examples on which the model has been trained. In addition, we have made this a
real-time system that captures ”snapshots” of the writer’s progress with
respect to the answer, allowing us to see how the student has arrived at
his/her final answer. This additional information can allow us to analyse
where students err, allowing us to better understand the thought-process of
students as they take tests and possibly unearth previously unknown trends.

##  I INTRODUCTION

Automated answer scoring is the use of specialized computer programs to assign
grades to answers written in an educational setting. It is a method of
educational assessment and an application of natural language processing.
Attempts to build an automated grading system dated back to 1966 when Ellis B.
Page proved on The Phi Delta Kappan that a computer could do as well as a
single human judge [1]. Since then, much effort has been put into building the
perfect system. Intelligent Essay Assessor (IEA), developed by Peter Foltz and
Thomas Landauer, was first used to score essays for large undergraduate
courses in 1994 [2]. The automated reader developed by the Educational Testing
Service, e-Rater, used hundreds of manually defined features. It was trained
on 64 different prompts and more than 25,000 answers on a 6-point scale from 1
to 6. Evaluated on the quadratic weighted kappa calculated between the
automated scores for the answers and the resolved score for human raters on
each set of essays, e-rater could only achieve a kappa score below 0.5. [3]
In 2012, the Hewlett Foundation sponsored a competition on Kaggle called the
Automated Student Assessment Prize (ASAP). The competition also used quadratic
weighted kappa to measure the similarity between the human scores and the
automated scores. 154 participants attempted to predict the essay score. The
winning team got a kappa score of 0.81407. [4] Later, a team at Carnegie
Mellon University built a model using dense and sparse features, trained on
the same dataset to achieve the kappa score of 0.833. [5]
Our approach will be different from the previous attempts to solve the
problem, using deep learning, as opposed to domain oriented predefined
features. The best way to improve one’s own writing skill is to write an
answer, receive feedback from your instructor, and based on the feedback,
revise your answer. Repetition of this process as often as possible is advised
up until an optimal score is reached. However, the problem arises in the fact
that continuously evaluating these essays requires a lot of time and puts an
enormous load on the classroom teacher who is often found in a position where
she must provide feedback to a possible 30 essays or more each time a topic is
assigned. As a result, teachers are not able to give writing assignments as
often as they would wish. Keeping this in mind, researchers have sought to
develop applications that automate essay scoring and evaluation, with the goal
of allowing teachers more flexibility in designing a course so as to maximize
the learning of the student without having to compromise due to a large work
load. Automated scoring of answers is not meant to serve as a replacement to a
human evaluator, but, merely as an aid to the evaluator.

##  II RELATED WORK

###  II-A Recurrent Neural Networks for Language Understanding [6]

Shows the effectiveness of Recurrent Neural Networks for language
understanding. The core of their project was to take words as inputs in a
standard RNN-LM and to predict slot labels instead of the traditional method
of words on the output side. The RNN-LM represents each word as a high-
dimensional real-valued vector in such a way that similar words tend to be
closer together and relationships between words are preserved. They
demonstrated the use of the Elman Architecture and using the ATIS Dataset,
they were able to show that the RNN outperformed the Conditional Random Fields
(CRF) approaches by a large margin as well as previous neural-net approaches.

###  II-B Automated Essay Scoring Using Bayes’ Theorem [7]

Bayesian Models for text classification from the information science field
were extended and applied to student produced essays, with each model
calibrated to 462 essays with no score points. They manipulated some
variables, such as trimming, stemming and the use of stopwords to imporve
their accuracy. They were able to achieve an accuracy of only 80%. The major
drawbacks of this paper are that it uses a very simple Bayesian Model for the
classification.

###  II-C Cosine similarity to determine similarity measure: Study case in
online essay assessment [8]

Implementation of the weighting of Term Frequency - Inverse Document Frequency
(TF-IDF) method and Cosine Similarity with the measuring degree concept of
similarity terms in a document. Tests have been carried out on a number of
Indonesian text-based documents that have gone through the stage of pre-
processing for data extraction purposes. This process results is in a ranking
of the document weight that have closeness match level with expert’s document.

##  III PROBLEM STATEMENT

The aim is to implement a system that evaluates a student’s answer $A$ and
recommends a score $s(A)$, accurately as well as efficiently, such that;
rather than just evaluating the end-product $A_{f}$, consider a sequence of
intermediate answers $A_{1}$, $A_{2}$, $...$, $A_{n}$ (where $A_{n}$ =
$A_{f}$) captured at suitable times $t_{1}$, $t_{2}$, $...$, $t_{n}$ as the
student builds his/her answer.
By relating $s(A_{i})$ vs $t_{i}$ and plotting the graph, we will be able to
gain insight on the thought process and steps taken by which the student
arrived at his/her answer. For instance, two students may end up with the same
score but one of the students may have gradually built up to this score,
whereas the other may have actually written an answer that earned a higher
score before making modifications that lowered the score.

###  III-A The Dataset

The data set we are using for this project is The Hewlett Foundation:
Automated Essay Scoring [4].

About the data set:

$\bullet$

There are 8 different sets of essays, each generated from a single prompt.

$\bullet$

Selected essays range from an average length of 150 to 550 words per response.

$\bullet$

Each essay was hand-graded a human evaluator, which means we are using real-
world data to produce real-world results.

$\bullet$

There are 10686 samples in total. We used a 85:15 data split for training and
testing, as well as validation.

$\bullet$

Each set has a different grading scale.

###  III-B Evaluation Metric [16]

In recent years, the role of big data analytics has exponentially grown, and
is now slowly making its way into the education industry. There is a growing
need to tap into the massive reserves of education data emanating from both
private as well as government institutions. Several attempts are being made in
this sphere in order to improve the quality of education being provided to
students while simultaneously improving teaching techniques to reduce the
burden put on the teacher. While the implications of such collaboration have
in fact been carried out before, automated scoring of answers has been
explored to a rather limited extent.

The biggest obstacle to choosing constructed-response assessments over the
traditional multiple-choice assessments is the large cost, and effort that
comes with their scoring, and this is precisely the issue that this project
aims to solve.The project is to accept raw-input from the student in the form
of their answer, pre-process these answers if required, and automatically
score the answers. Evaluation of these answers will be based purely on the
previous examples on which the model has been trained. In addition, we have
made this a real-time system that captures ”snapshots” of the writer’s
progress with respect to the answer, allowing us to see how the student has
arrived at his/her final answer. This additional information can allow us to
analyse where students err, allowing us to better understand the thought-
process of students as they take tests and possibly unearth previously unknown
trends.

κ= 1 - ∑ijWijOij∑ijWijEij

$\bullet$

The Weight matrix $W$ is calculated based on the difference between raters
scores.

$\bullet$

$O_{ij}$ corresponds to the number of answers that received a rating $i$ by
Rater A and Rating $j$ by Rater B.

$\bullet$

$E$ is the $nxn$ histogram matrix of expected ratings.

The evaluation metric used for the project is quadratic weighted kappa [10].
The Quadratic Weighted Kappa (QWK) metric typically varies from 0 (only random
agreement between raters) to 1 (complete agreement between raters), and
doesn’t penalize for an incorrect score match, but rather, it takes into
account the error in the scores, with respect to the range of scores.
Therefore, it serves as a much better evaluation metric as compared to
accuracy, since in our problem, we’re primarily comparing inter-rater
agreement.

##  IV APPROACH

The primary objective is to build a model that accepts an answer as input and
automatically outputs the score of that answer, in an efficient manner that
can be done in real-time as the student writes the answer. The score would be
evaluated at the end of sentences, words or ever a few systems that only
perform the analysis and evaluate the final draft in a answer.
Coming to the pre-processing don. What we first did, was try to build a Bag of
Words Model and examine the results it would give us, but they were not very
promising. Based on our findings, we decided to move forward and remove all
stop-words from each of the essays being passed to our Model. We also
implemented a rudimentary Spell Check on the essays based on the Peter Norvig
Spell Check Algorithm and stemming, however, the results we achieved were
substandard possibly owing to the fact that human evaluators would consider a
misspelled word as an important factor (reducing overall marks of the answer)
and we decided against using them continuing forward with our research. The
interesting bit of the pre-processing comes, where we decided to implement a
word-embedding model. As mentioned above we had two possible choices, GloVe
[12] and Word2Vec [16]. Based on some research into the capabilities of both
models we found that the Word2Vec would possibly be a better model, however,
just to confirm for ourselves, we decide to implement both, working however,
on the assumption that the Word2Vec Model would outperform the GolVe Model.
The role that Word2Vec plays in our implementation is a massive one. Developed
by team of researchers at Google, it is an extremely useful group of related
models to produce word vectors of your input. The Word2Vec Model is trained on
the input specified, and is trained in such a way that words that have a
similar meaning have a smaller vector distance between them. This
representation of words is very helpful going forward with our predictions
primarily because we are able to train our Model based on the type of input,
rather than predefined meanings of the word, which was consequently one of the
reasons that we chose this model over GloVe, thereby ensuring that words that
have a similar meaning in that context or tend to appear together frequently
have smaller distances than words that are similar by their intrinsic meaning.

TABLE I: Random Forest v/s SVM RBF | Word2Vec |
---|---|---
Model | QWK | Dim
Random Forests | 0.9535 | 300
Random Forests | 0.9563 | 200
Random Forests | 0.9508 | 100
SVM Radial Basis Function | 0.9619 | 300
SVM Radial Basis Function | 0.9568 | 200
SVM Radial Basis Function | 0.9509 | 100

Once we obtained the vector representations of our input, we used Vector
Averaging to get a specific Vector Value for each essay in the dataset, which
we could then use to make better predictions. In our approach we used four
Models; a Random Forest Classifier, a Radial Basis Support Vector Machine, a
Deep Neural Network and a Short Long Term Memory Unit which is a form of a
Recurrent Neural Network. The Random Forest Classifier [17] and the Support
Vector Machine [18] returned very good results. Looking into Random Forest
Classifiers, they are an example of an ensemble learner built on decision
trees, is a direct consequence of the fact that by maximum voting from a panel
of independent judges, we get the final prediction better than the best judge.
Our Deep Neural Network is a simple two-layered feed forward neural network,
where we use the ’Rectified Linear Units’ as the activation function for the
first layer and a ’Softmax’ activation function [19] for the output layer to
squash the output vector. We tried multiple activation functions, but verified
that these activation functions give us an optimum result and the loss
function that we have used is ’sparse categorical crossentropy’ since our
problem deals with more than two classes. While building the Word2Vec Model we
achieved the best results when the dimension for each word vector was set to
300. Once trained, the Neural Network architecture would be used to evaluate
’snapshots’ of the students’ intermediate responses which will give us a graph
of the student’s score vs time intervals, which could answer many interesting
questions with respect to the student’s performance.

TABLE II: LSTM v/s DNN v/s Combined Model | Word2Vec | GloVe |  |
---|---|---|---|---
Model | QWK | QWK | Dim |
Long Short Term Memory | 0.9653 | 0.9566 | 300 |
Long Short Term Memory | 0.9642 | 0.9623 | 200 |
Long Short Term Memory | 0.9512 | 0.9433 | 100 |
Deep Neural Network | 0.9643 | 0.9555 | 300 |
Deep Neural Network | 0.9576 | 0.955 | 200 |
Deep Neural Network | 0.9611 | 0.9442 | 100 |
Combined Model | 0.9721 | - | 100 |

Recurrent networks take as their input not just the current input example they
see, but also what they have perceived previously in time. The decision a
recurrent net reached at time step t-1 affects the decision it will reach one
moment later at time step t. Hence, RNNs have two sources of input, the
present and the recent past, which combine to determine how they respond to
new data. This quality of an RNN is what allows it to produce robust models
when it comes to textual input, because of its sequential quality. The
sequential information is preserved in the RNN’s hidden state, which manages
to span many time steps as it cascades forward to affect the processing of
each new example. It is finding correlations between events separated by many
moments, and these correlations are called “long-term dependencies”, because
an event downstream in time depends upon, and is a function of, one or more
events that came before.
In the mid-90s, a variation of recurrent net, called Long Short Term Memory
units (LSTM) was proposed as a solution to the vanishing gradient problem
[24]. LSTMs help preserve the error that can be back-propagated through time
and layers. By maintaining a more constant error, they allow RNNs to continue
to learn over many time steps (over 1000), thereby opening a channel to link
causes and effects remotely. This feature of an LSTM, aids the RNN in the
model we have built to allow for better predictions of each answer.
Now, considering that we are trying to predict a grade for a student, a simple
accuracy test where we compare the number of correct predictions of our model,
would not be a very good metric, therefore, as specified in the source of the
dataset, have used Quadratic Weighted Kappa as our Evaluation Metric.
Quadratic Weighted Kappa allows us to calculate the agreement between
evaluator, which in our case would be the agreement of the score predicted by
our models, against the score predicted by an evaluator whose scores we are
using to train the model. This turns out to be a model far more similar to a
real-world scenario where agreement between raters is considered, taking into
account the weighted errors.
Finally, we used a method similar to Weighted Average. We obtained the
predictions from each Model that we built, namely the LSTM Mode, the Deep
Neural Network Mode, the SVM Radial Basis Function Model and the Random
Forrest Classifier Model, along with their respective Quadratic Weighed Kappa
values, and applying the formula:

κ_f = κ1*S1\+ κ2*S2\+ κ3*S3\+ κ4*S4∑i=14κi

we were able to obtain better predictions for each essay that resulted in an
overall higher Quadratic Weighted Kappa Score.

###  IV-A General Classifier Results

Table 2 above shows us the QWK of the Random Forrest Classifier and the SVM
using a Radial Basis Function. The high QWK was unexpected, however we
attribute it towards the pre-processing done as well as the use of Word
Vectors.

##  V EXPERIMENTS AND RESULTS

Based on the above approaches, we ran multiple tests to figure out which Model
would give the best results. We documented all the results and we now provide
justifications for each of them, along with a couple other pre-processing
techniques that we used, and how they fared as well. Earlier, we had justified
our selection of Word2Vec over the pre-trained GloVe Vectors, but as a
confirmation, we ran tests on both models just to take a look at how they
performed under the same conditions.

###  V-A Effect of Pre-processing

We used a couple of different pre-processing techniques which affected our
results, so we feel it is appropriate to discuss these before our final
Results

####  V-A1 Stemming [20]

Stemming is the process of reducing inflected words to their word stem, base
or root form. When we applied this method we incurred a drop in the QWK across
all models. The reason behind this drop lies in the importance of the way in
which an answer is written by a student. A human scorer would evaluate
students who have better grammar and use of words higher than those who did
not. Hence, when we tested the model after stemming the training and testing
sets, we incurred a considerable drop in accuracy.

####  V-A2 Spell-Correction [11]

To implement Spell-Correction we used Peter Norvig’s Spell Check Algorithm
[11], and this method as well, similarly to Stemming, resulted in a drop of
accuracy. The reason for this is identical to the one listed above, and hence,
we deiced not to use either of these methods moving forward with our testing.

Figure 1: Comparison of the different models

###  V-B Neural Network Results

Table 1 above, has a comparison of how the LSTM Model and Deep Neural Network
Model worked with respect to both Word Embedding Models. Based on the results,
the following two observations can be made:

####  V-B1 Word Embedding Model

You can see that in general the Word2Vec Model outperforms the GloVe Model.
This is because Word2Vec was trained on our own dataset whereas the GloVe
Model used pre-trained Word Vectors. This phenomenon is best explained with an
example. Let’s say we have two words, threads and string, in terms of
Computers (which was the topic of most of the answers in our dataset) they
have a completely different meaning to what they do in day-to-day life. Hence,
while predicting a score, the Word2Vec Model allows us to use the relevant
word in the relevant situation which GloVe does not, leading to a slightly
lower Quadratic Weighted Kappa.

####  V-B2 Prediction Model

We see that the LSTM outperforms the DNN, this increase in QWK can be
attributed towards the fact that an LSTM is basically a Recurrent Neural
Network, this means that it works best with sequential data, and in our case,
text is a perfect example of sequential data. Additionally, an LSTM has a
memory unit which allows it to take into account a larger context, thereby
giving slightly better results. However we see that the Combined Model yields
the best overall results, due to the fact that it takes into account the
prediction made by each Model and gives a weighted average that tends to be
the best fit for the Answer in question.

Figure 2: Errors of the Combined Model using Word2Vec

###  V-C Graphs

The visualizations seen is a general error graph of the predictions obtained
using the Combined Model.
Finally you can see our results compared to the previous best results of the
competitors who participated in Kaggle’s competition.
We also added a graph of QWK v/s the Model used which is a comparison of all
our Models compared to previous results. As illustrated in the graph, our
Model has exceeded previous prediction Model by a large margin.

###  V-D Optimizers

With regard to optimizers, we used a numerous predefined optimizers such as
$StochasticGradientDescent$ [21], $Adam$ [22], $Adagrad$, $AdaDelta$ [23] etc.
But $RMSProp$ yielded the best results, with a Quadratic Weighted Kappa of
about 0.974

###  V-E Real-Time System

We implemented our Real-Time system using a web-server designed on python’s
flask web framework. We grade the scores of students as they type at intervals
of 0.66 seconds based on studies that show this to be the typing speed per
word of the average human. The Real-time system predicts the score of
”snapshots” of a student’s answer using the combined model in about 0.3
seconds (running on a Intel i7-7820 Processor), making it an interactive
system.

##  VI Conclusion

The inspiration behind this project was to make the job of a teacher easier.
The crux of a formal education is written assignments, and teachers often are
not able to give as many as they would like. Using a system built around
Automated Scoring, we could move into an era where teachers can Model courses
however they wish without being restricted by their workload.
We also find that this project could be helpful in getting to know how a
student answers particular questions, and in what way they first approach a
particular problem. This would give us a huge amount of data with respect to
how students think, this giving us a better idea on how to teach certain
concepts to students, in ways that they best understand.

## References

  * [1] Page, Ellis B. ”The imminence of… grading essays by computer.” The Phi Delta Kappan 47.5 (1966): 238-243..
  * [2] Foltz, Peter W., et al. ”Implementation and applications of the Intelligent Essay Assessor.” Handbookof automated essay evaluation (2013): 68-88.
  * [3] Attali, Yigal, and Jill Burstein. ”Automated essay scoring with e-rater V. 2.” The Journal of Technology, Learning and Assessment 4.3 (2006).
  * [4] Kaggle. ”Develop an automated scoring algorithm for student-written essays.” (2012). https://www.kaggle.com/c/asap-aes
  * [5] Robertson, Stephen. ”Understanding inverse document frequency: on theoretical arguments for IDF.” Journal of documentation 60.5 (2004): 503-520.
  * [6] Yao, Kaisheng, Geoffrey Zweig, Mei-Yuh Hwang, Yangyang Shi, and Dong Yu. ”Recurrent neural networks for language understanding.” In Interspeech, pp. 2524-2528. 2013.
  * [7] Rudner, Lawrence M., and Tahung Liang. ”Automated essay scoring using Bayes’ theorem.” The Journal of Technology, Learning and Assessment 1, no. 2 (2002).
  * [8] Lahitani, Alfirna Rizqi, Adhistya Erna Permanasari, and Noor Akhmad Setiawan. ”Cosine similarity to determine similarity measure: Study case in online essay assessment.” In Cyber and IT Service Management, International Conference on, pp. 1-6. IEEE, 2016.
  * [9] Foltz, Peter W., Darrell Laham, and Thomas K. Landauer. ”Automated essay scoring: Applications to educational technology.” In World Conference on Educational Multimedia, Hypermedia and Telecommunications, vol. 1, pp. 939-944. 1999.
  * [10] Cohen, Jacob. ”Weighted kappa: Nominal scale agreement provision for scaled disagreement or partial credit.” Psychological bulletin 70, no. 4 (1968): 213.
  * [11] Norvig, Peter. ”Natural language corpus data.” Beautiful Data (2009): 219-242.
  * [12] Pennington, Jeffrey, Richard Socher, and Christopher Manning. ”Glove: Global vectors for word representation.” In Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP), pp. 1532-1543. 2014.
  * [13] Ramos, Juan. ”Using tf-idf to determine word relevance in document queries.” In Proceedings of the first instructional conference on machine learning, vol. 242, pp. 133-142. 2003.
  * [14] Sepp Hochreiter and Jürgen Schmidhuber, Long Short-Term Memory, Neural Computation Volume 9, Issue 8, November 15, 1997 p.1735-1780
  * [15] Mikolov, Tomas; et al. ”Efficient Estimation of Word Representations in Vector Space”. https://arxiv.org/abs/1301.3781
  * [16] Mary L. McHugh, Inter-Rater Reliability: The Kappa Statistic.
  * [17] Breiman, L. Machine Learning (2001) 45: 5. https://doi.org/10.1023/A:101093340432
  * [18] M.A. Hearst ; S.T. Dumais ; E. Osuna ; J. Platt ; B. Scholkopf, Support Vector Machines Published in: IEEE Intelligent Systems and their Applications (Volume: 13, Issue: 4, July-Aug. 1998)
  * [19] R. A. Dunne and N. A. Campbelly, ”On The Pairing Of The Softmax Activation And Cross Entropy Penalty Functions And The Derivation Of The Softmax Activation Function” 8th Aust.Conf.on Neural Networks, Melb.1997
  * [20] Julie Beth Lovins, Development of a Stemming Algorithm, [Mechanical Translation and Computational Linguistics, vol.11, nos.1 and 2, March and June 1968]
  * [21] Antoine Bordes, Léon Bottou, Patrick Gallinari, SGD-QN: Careful Quasi-Newton Stochastic Gradient Descent Journal of Machine Learning Research 10 (2009) 1737-1754 Submitted 1/09; Revised 4/09; Published 7/09
  * [22] Diederik P. Kingma, Jimmy Ba, Adam: A Method for Stochastic Optimization
  * [23] Matthew D. Zeiler, ADADELTA: An Adaptive Learning Rate Method
  * [24] Y. Bengio, P. Simard, P. Frasconi: Learning long-term dependencies with gradient descent is difficult


Department of Computer Science, University of Helsinki, Finland and
diediaz@dcc.uchile.clhttps://orcid.org/0000-0002-1825-0097Academy of Finland
Grant 323233 CeBiB – Centre For Biotechnology and Bioengineering, Department
of Computer Science, University of Chile,
<EMAIL_ADDRESS>Basal Funds FB0001 and Fondecyt Grant
1-200038, Chile Diego Díaz and Gonzalo Navarro [500]Theory of computation Data
compression John Q. Open and Joan R. Access 2 42nd Conference on Very
Important Topics (CVIT 2016) CVIT 2016 CVIT 2016 December 24–27, 2016 Little
Whinging, United Kingdom 42 23

# Efficient Construction of the BWT for Repetitive Text Using String
Compression

Diego Díaz-Domínguez  Gonzalo Navarro

###### Abstract

We present a new semi-external algorithm that builds the Burrows–Wheeler
transform variant of Bauer et al. (a.k.a., BCR BWT) in linear expected time.
Our method uses compression techniques to reduce the computational costs when
the input is massive and repetitive. Concretely, we build on induced suffix
sorting (ISS) and resort to run-length and grammar compression to maintain our
intermediate results in compact form. Our compression format not only saves
space, but it also speeds up the required computations. Our experiments show
important savings in both space and computation time when the text is
repetitive. On average, we are 3.7x faster than the baseline compressed
approach, while maintaining a similar memory consumption. These results make
our method stand out as the only one (to our knowledge) that can build the BCR
BWT of a collection of 25 human genomes (75 GB) in about 7.3 hours, and using
only 27 GB of working memory.

###### keywords:

BWT, string compression, repetitive text

###### category:

##  1 Introduction

The Burrows–Wheeler transform (BWT) [6] is a reversible string transformation
that reorders the symbols of a text $T$ according the lexicographical ranks of
its suffixes. The features of this transform have turned it into a key
component for text compression and indexing [32, 25]. In addition to being
reversible, the reordering produced by the BWT reduces the number of equal-
symbol runs in $T$, thus improving the compresibility. On the other hand, its
combinatorial properties [10] enable the creation of self-indexes [9, 28] that
support pattern matching in time proportional to the pattern length. Popular
bionformatic tools [18, 21] rely on the BWT to process data, as collections in
this discipline are typically massive and repetitive, and the patterns to
search for are short.

There are several algorithms in literature that produce the BWT in linear time
[34, 1, 23, 8, 3]. Nevertheless, the computational resources their
implementations require when the input is large are still too high for
practical purposes. This problem is particularly evident in Genomics
applications, where the amount of data is growing at an astronomical rate
[35].

Although genomic collections are becoming more and more massive, the effective
information they contain remains low compared to their sizes [27]. A promising
solution to deal with this kind of data is then to design BWT algorithms that
scale with the amount of information in the collection, not with its size.

Motivated by these ideas, some authors have developed new BWT algorithms that
exploit text repetitions to reduce the computational requirements [14, 5, 13,
15, 4]. Their approach consists of extracting a set of representative strings
from the text, perform calculations on them, and then extrapolate the results
to the copies of those strings. For instance, the methods of Boucher et al.
[5, 4] based on prefix-free parsing (PFP) use Karp–Rabin fingerprints [12] to
create a dictionary of prefix-free phrases from $T$. Then, they create a parse
by replacing the phrases in $T$ with metasymbols, and finally construct the
BWT using the dictionary and the parse. Similarly, Kempa et al. [14] consider
a subset of positions in $T$ that they call a string synchronizing set, from
which they compute a partial BWT they then extrapolate to the whole text.

Although these repetition-aware techniques are promising, some of them are at
a theoretical stage [14, 13, 15], while the rest [5, 4] have been empirically
tested only under controlled settings, and their results depend on parameters
that are not simple to tune. Thus, it is difficult to assess their performance
under real circumstances.

Recently, Nunes et al. [31] proposed a method called GCIS that adapts the
concept of _induced suffix sorting_ (ISS) for compression. Their ideas are
closely related to the linear-time BWT algorithm of Okanohara et al. [34].
Briefly, Okanohara et al. cut the text into phrases using ISS, assign symbols
to the phrases, and then replace the phrases with their symbols. They apply
this procedure recursively until all the symbols in the text are different.
Then, when they go back from the recursions, they induce an intermediate
$BWT^{i}$ for the text of every recursion $i$ using the previous $BWT^{i+1}$.
On the other hand, GCIS stores the dictionaries that ISS generates in the
recursions in a context-free grammar. The connection between these two methods
is that GCIS captures in the grammar precisely the information that Okanohara
et al. use to compute the BWT. Additionally, Díaz-Domínguez et al. [7]
recently demonstrated that ISS-based compressors such as GCIS require much
less computational resources than state-of-the-art methods like RePair [19] to
encode the data, while maintaining high compression ratios. The simple
construction of ISS makes it an attractive alternative to process high volumes
of text. In particular, combining the ideas of Okanohara et al. with ISS-based
compression is a promising alternative for computing big BWTs.

Our contribution. Induced suffix sorting (ISS) [17] has proved useful for
compression [31, 7] and for constructing the BWT [34]. In this work, we show
that compression can be incorporated in the internal stages of the BWT
computation in a way that saves both working space and time. Okanohara et al.
[34] use ISS to construct the BWT in recursive stages of parsing (which cuts
the text into metasymbols) and partial construction of the BWT of the
metasymbols; the final BWT is obtained when returning from the recursion. We
use a technique similar to grammar compression to store the dictionaries of
metasymbols, and run-length compression for the partial BWTs. This approach is
shown not only to save the space required for those intermediate results, but
importantly, the format we choose actually speeds up the computation of the
final BWT as we return from the recursion, because the factorizations that
help save space also save redundant computations. Unlike Okanohara et al., we
receive as input a string collection and output its BCR BWT [1], a variant for
string collections. The reason is that massive datasets usually contain
multiple strings, in which case the BCR BWT variant is simpler to construct.
Our experiments show that, when the input is a collection of human genomes (a
repetitive dataset), our implementation requires 3.7x less computation time
than ropeBWT2 [21], an efficient implementation of the BCR BWT algorithm.
Additionally, we use 7.6x less working memory than pfp-ebwt [4], a recent
method that uses an strategy similar to ours. Under not so repetitive
scenarios, our performance is competitive with BCR_LCP_GSA [1] or gsufsort
[8].

##  2 Related Concepts

###  2.1 The Burrows–Wheeler Transform

Consider a string $T[1,n-1]$ over alphabet $\Sigma[2,\sigma]$, and the
sentinel symbol $\Sigma[1]=\texttt{\$}$, which we insert at $T[n]$. The
_suffix array_ [26] of $T$ is a permutation $S\\!A[1,n]$ that enumerates the
suffixes $T[i,n]$ of $T$ in increasing lexicographic order,
$T[S\\!A[i],n]<T[S\\!A[i+1],n]$, for $i\in[1,n-1]$.

The _Burrows–Wheeler transform_ (BWT) [6] is a reversible string
transformation that stores in $BWT[i]$ the symbol that precedes the $ith$
suffix of $T$ in lexicographical order, i.e., $BWT[i]=T[S\\!A[i]-1]$ (assuming
$T[0]=T[n]=\texttt{\$}$).

The mechanism to revert the transformation is the so-called $\mathsf{LF}$
mapping. Given an input position $BWT[j]$ that maps a symbol $T[i]$,
$\mathsf{LF}(j)=j^{\prime}$ returns the index $j^{\prime}$ such that
$BWT[j^{\prime}]=T[i-1]$ maps the preceding symbol of $T[i]$. Thus, spelling
$T$ reduces to continuously applying LF from $BWT[1]$, the symbol to the left
of $T[n]=\texttt{\$}$, until reaching $BWT[j]=\texttt{\$}$.

The BCR BWT [1] is a variant of the original BWT that encodes a string
collection $\mathcal{T}=\\{T_{1},T_{2},\ldots,T_{k}\\}$ instead of a single
string $T$. Briefly, if two (or more) symbols $a=T_{x}[k]$ and
$b=T_{y}[k^{\prime}]$, from different strings $T_{x},T_{y}\in\mathcal{T}$, are
preceded by identical suffixes $T_{x}[k+1..]=T_{y}[k^{\prime}+1..]$, the order
of $a$ and $b$ in $BWT$ is the same as the relative order of $T_{x}$ and
$T_{y}$ in $\mathcal{T}$. The BCR BWT also appends sentinel symbols $ to the
strings of $\mathcal{T}$ to detect their boundaries in the BWT. A position
$BWT[j]=\texttt{\$}$ represents the start of a string $T_{u}$, and
$BWT[\mathsf{LF}(j)]$ maps the end of a string $T_{u^{\prime}}\in\mathcal{T}$
that is not necessarily $T_{u}$.

###  2.2 Grammar and Run-length Compression

_Grammar compression_ [16] consists of encoding a text $T$ as a small context-
free grammar $\mathcal{G}$ that only produces $T$. Formally, a grammar is a
tuple $(V,\Sigma,\mathcal{R},\mathqhv{S})$, where $V$ is the set of
nonterminals, $\Sigma$ is the set of terminals, $\mathcal{R}$ is the set of
replacement rules and $\mathqhv{S}\in V$ is the start symbol. The right-hand
side of $\mathqhv{S}\rightarrow C\in\mathcal{R}$ is referred to as the
compressed form of $T$. The size of $\mathcal{G}$ is usually measured in terms
of the number of rules, the sum of the lengths of the right-hand sides of
$\mathcal{R}$, and the length of the compressed string.

_Run-length compression_ encodes the equal-symbol runs of maximal length in
$T$ as pairs. More specifically, $T$ becomes a sequence
$(a_{1},l_{1}),(a_{2},l_{2}),\ldots,(a_{n^{\prime}},l_{n^{\prime}})$ of
$n^{\prime}\leq n$ pairs, where every $(a_{i},l_{i})$, with
$i\in[1,n^{\prime}]$, stores the symbol $a_{i}\in\Sigma$ of the $ith$ run and
its length $l_{i}\geq 1$. For instance, let $T[i,j]=aaaa$ be a substring with
four consecutive copies of $a$, where $T[i-1]\neq a$ and $T[j+1]\neq a$. Then
$T[i,j]$ compresses to $(a,4)$.

###  2.3 Induced Suffix Sorting

_Induced suffix sorting_ (ISS) [17] computes the lexicographical ranks of a
subset of suffixes in $T$ and then uses the result to induce the order of the
rest. This method is the underlying procedure in several algorithms that build
the suffix array [30, 29, 22] and the BWT [34, 5] in linear time. The ISS idea
introduced by the suffix array algorithm SA-IS of Nong et al. [30] is of
interest to this work. The authors give the following definitions:

######  Definition 2.1.

A symbol $T[i]$ is called L-type if $T[i]>T[i+1]$ or if $T[i]=T[i+1]$ and
$T[i+1]$ also L-type. On the other hand, $T[i]$ is said to be S-type if
$T[i]<T[i+1]$ or if $T[i]=T[i+1]$ and $T[i+1]$ is also S-type. By definition,
symbol $T[n]$, the one with the sentinel, is S-type.

######  Definition 2.2.

A symbol $T[i]$, with $i\in[1,n]$, is called leftmost S-type, or LMS-type, if
$T[i]$ is S-type and $T[i-1]$ is L-type.

######  Definition 2.3.

An LMS substring is (i) a substring $T[i,j]$ with both $T[i]$ and $T[j]$ being
LMS symbols, and there is no other LMS symbol in the substring, for $i\neq j$;
or (ii) the sentinel itself.

SA-IS is a recursive method. In every recursion $i$, it initializes an empty
suffix array $A^{i}$ for the input text $T^{i}$ (i=1). Then, it scans $T^{i}$
from right to left to classify the symbols as L-type, S-type or LMS-type. As
it moves through the text, the algorithm records the text positions of the LMS
substrings in $A^{i}$. More specifically, if $T^{i}[j]=a$ is the first symbol
of an LMS substring, it inserts $j$ in the right-most empty position in the
bucket $a$ of $A^{i}$. After scanning $T^{i}$, SA-IS sorts the LMS substrings
in $A^{i}$ using ISS. This procedure only requires two linear scans of $A^{i}$
(we refer the reader to Nong et al. [30] for further detail).

ISS sorts the LMS substrings in a way that is slightly different from
lexicographic ordering, we refer to it as $\prec_{LMS}$ ordering. In
particular, if an LMS substring $T^{i}[a,b]$ is a prefix of another LMS
substring $T^{i}[a^{\prime},b^{\prime}]$, then $T^{i}[a,b]$ gets higher order.
However, the higher rank of $T^{i}[a,b]$ implies that the suffix $T^{i}[a..]$
is lexicographically greater than the suffix $T^{i}[a^{\prime}..]$. The cause
of this property is explained in Section 2 of Ko and Aluru [17].

The idea now is to use the sorted LMS substrings to induce the order of the
suffixes in $T^{i}$ that are not prefixed by LMS substrings. Still, LMS
substring with the same sequence are still unsorted in $A^{i}$. Nong et al.
solve this problem by creating a new string $T^{i+1}$ in which they replace
the distinct LMS substrings with their orders in $A^{i}$, and use $T^{i+1}$ as
input for another recursive call $i+1$. The base case for the recursion is
when all the suffixes in $A^{i}$ are prefixed by different symbols, in which
case they return $A^{i}$ without further processing.

When the $(i+1)th$ recursive call ends, the suffixes of $T^{i}$ prefixed by
the same LMS substrings are completely sorted in $A^{i+1}$, so SA-IS proceeds
to complete $A^{i}$. For doing so, it resets $A^{i}$, inserts the LMS
substrings arranged as they respective symbols appear in $A^{i+1}$, and
performs ISS again to reorder the unsorted suffixes of $T^{i}$. Once it
finishes, it passes $A^{i}$ to the previous recursion $i-1$. The final array
$A^{1}$ is the suffix array for $T$.

##  3 Methods

###  3.1 Definitions

We consider a collection $\mathcal{T}=\\{T_{1},T_{2},\ldots,T_{k}\\}$ of $k$
strings over the alphabet $\Sigma[2,\sigma]$. The input for our algorithm is
thus the sequence $T=T_{1}\texttt{\$}T_{2}\ldots T_{k}\texttt{\$}$ of total
length $n=|T|$ that represents the concatenation of $\mathcal{T}$. The symbol
$ is a sentinel that we use as a boundary between consecutive strings in $T$.
We map $\texttt{\$}=\Sigma[1]$ to the smallest symbol in the alphabet.

Let $\mathcal{D}=\\{D_{1},D_{2},\ldots,D_{y}\\}$ be a string set that is not
suffix-free. A suffix $D_{j}[u..]$, with $D_{j}\in\mathcal{D}$, is _proper_ if
$1<u\leq|D_{j}|$. Additionally, we consider a suffix $D_{j}[u..]$, with
$D_{j}\in\mathcal{D}$, to be _left-maximal_ if there is at least one other
suffix $D_{j^{\prime}}[u^{\prime}..]$, with $D_{j^{\prime}}\in\mathcal{D}$,
such that (i) $j\neq j^{\prime}$, (ii)
$D_{j^{\prime}}[u^{\prime}..]=D_{j}[u..]$, and (iii) both
$D_{j^{\prime}}[u^{\prime}..]$ and $D_{j}[u..]$ are proper suffixes with
$D_{j^{\prime}}[u^{\prime}-1]\neq D_{j}[u-1]$ or one of them is not a proper
suffix.

###  3.2 Overview of Our Algorithm

We call our algorithm for computing the BCR BWT of $\mathcal{T}$
$\mathsf{grlBWT}$. This method relies on the ideas developed by Nong et al. in
the SA-IS algorithm (Section 2.3), but includes elements of grammar and run-
length compression (Section 2.2). These new features reduce the space usage of
the temporal data that $\mathsf{grlBWT}$ maintains in memory, thus decreasing
both working memory and computing time. We now give a brief overview of our
approach.

Our method works in two phases: the _parsing_ phase and the _induction_ phase.
The parsing phase is similar to the recursive steps of SA-IS. In every
iteration $i$ (or parsing round), we first scan the input string $T^{i}$
($T^{1}=T$) to build a dictionary $D^{i}$ with the phrases that occur as LMS
substrings. We also record the frequency of every phrase, i.e., the number of
times it occurs as an LMS substring in $T^{i}$. Subsequently, we use the
phrases in $D^{i}$ and their frequencies to construct a preliminary BWT for
$T^{i}$ ($pBWT^{i}$), which we complete in the induction phase. We say
$pBWT^{i}$ is partial because it has empty spaces we can not fill just with
the information in $D^{i}$. To make the completion more efficient during the
induction phase, we encode $pBWT^{i}$ using run-length compression and $D^{i}$
using a technique similar to grammar compression. Finally, we store $pBWT^{i}$
and $D^{i}$ on disk, and create a new text $T^{i+1}$ for the next parsing
round $i+1$. We construct $T^{i+1}$ by replacing the LMS substrings of $T^{i}$
with their associated symbols in $D^{i}$. The parsing phase finishes when no
new dictionary phrases can be extracted from the input text $T^{i}$ (see
Section 3.3).

Let $h$ be the number of iterations the parsing phase of $\mathsf{grlBWT}$
incurred with $T$. The induction phase starts by building the BWT for $T^{h}$.
After obtaining $BWT^{h}$, we start a new iterative process in which we
revisit the data we dumped to disk during the parsing phase in reverse order
(i.e., from round $h-1$ to round $1$). In every iteration $i$, the $BWT^{i+1}$
of $T^{i+1}$ is already computed, and we use it along with compressed version
of $D^{i}$ to induce the order of the symbols in the empty entries of
$pBWT^{i}$. Once we finish the induction, we compact $pBWT^{i}$ using run-
length encoding to create $BWT^{i}$. The final BCR BWT for $T$ is thus in
$BWT^{1}$.

###  3.3 The Parsing Phase

In this section, we explain the steps we perform during the $ith$ iteration of
the parsing phase of $\mathsf{grlBWT}$. Assume we receive as input a text
$T^{i}$ over the alphabet $\Sigma^{i}=[1,\sigma^{i}]$. We first initialize a
hash table $H^{i}$ that we will use to construct the dictionary $D^{i}$. The
keys in $H^{i}$ will be the phrases that occur as LMS substrings in $T^{i}$
while the values of $H^{i}$ will be the frequencies of the keys, i.e., the
number of times the keys occur as LMS substrings in $T^{i}$.

The basic idea to fill $H^{i}$ consists of scanning $T^{i}$ from right to left
to classify its symbols according the definitions of Section 2.3, and hash an
LMS subtring every time we reach an LMS-type symbol. This mechanism is almost
the same as the one described by Nong et al. [30] to detect the LMS substrings
(except for the hashing). However, we add an extra consideration. The
detection of LMS substrings is oblivious of the fact that $T^{1}$ encodes a
string collection rather than a single string. More precisely, in any parsing
round $i>1$, we could have an LMS substring of $T^{i}$ whose _expansion_
111Let $T^{i}[j,j^{\prime}]$ be a substring of $T^{i}$. We define the
_expansion_ of $T^{i}[j,j^{\prime}]$ as the string in $\Sigma^{1}$ we obtain
by recursively replacing every symbol $T^{i}[j]\in\Sigma^{i}$, for
$j\in[j,j^{\prime}]$, with its corresponding phrases in the dictionaries
$D^{i-1},D^{i-2},\ldots,D^{1}$. produces a substring of $T^{1}$ that covers
two or more strings of $\mathcal{T}$. These border phrases make the
computation of the BCR BWT a bit more difficult as we need to treat them
differently. We avoid this problem by maintaining a bit vector
$B^{i}[1,\sigma^{i}]$ that marks which symbols in $T^{i}$ expand to suffixes
of strings in $\mathcal{T}$. Thus, during the right-to-left scan of $T^{i}$,
each time we reach a position $T^{i}[j]$, such that $B^{i}[T^{i}[j]]=1$, we
truncate the active LMS substring. We record the phrase
$F=T^{i}[j+1,j^{\prime}]$ in $H^{i}$, where $T[j^{\prime}]$ is the last LMS-
type symbol we accessed, and start a new phrase from position $T^{i}[j]$.
Notice that the truncated strings are not LMS substrings by definition, but
they do not affect our algorithm (the reasons are explained in Definition 4
and Lemma 3 of Díaz-Domínguez et al. [7]). We receive $B^{i}$ as input along
with $T^{i}$ at the beginning of the parsing round $i$, and we compute the
next $B^{i+1}[1,\sigma^{i+1}]$ when we finish the round.

For practical reasons, we change the representation of $D^{i}$, encoded in
$H^{i}$ for the moment, to a more convenient data structure. First, we
concatenate all the keys of $D^{i}$ in one single vector $R^{i}$. We mark the
boundaries of consecutive phrases in $R^{i}$ with a bit vector $L^{i}$ in
which we set $L^{i}[j]=1$ if $R^{i}[j]$ is the first symbol of a phrase, and
set $L^{i}[j]=0$ otherwise. We also augment $L^{i}$ with a data structure that
supports $\mathsf{rank}_{1}$ queries [33] to map each symbol $D^{i}[j]$ to its
corresponding phrase. We store the values of $D^{i}$ in another vector
$N^{i}[1,|D^{i}|]$. We maintain the relative order so that the value
$N^{i}[o]$ maps the $oth$ phrase we inserted into $R^{i}$. For simplicity, we
will refer to the representation $(R^{i},L^{i},N^{i})$ just like $D^{i}$. We
still need $H^{i}$ to construct the parse $T^{i+1}$, so we do not discard it
but store it into disk.

The next step is to build $pBWT^{i}$ from $D^{i}$. For that purpose, we use
the following observations:

######  Lemma 3.1.

Let $X[1,x]$ and $Y[1,y]$ be two different strings over the alphabet
$\Sigma^{i}$, with lengths $x>1$ and $y>1$ (respectively). Assume both occur
as suffixes in one or more phrases of $D^{i}$. Let $\mathcal{X}$ be the list
of positions in $T^{i}$ where $X$ occurs as a suffix of an LMS substring. More
specifically, each $j\in\mathcal{X}$ is a position such that $T^{i}[j,j+x-1]$
is an occurrence of $X$ and $T^{i}[j-j^{\prime},j+x-1]$, with $j^{\prime}\geq
0$, is an LMS substring. Let us define a list $\mathcal{Y}$ equivalent to
$\mathcal{X}$, but for $Y$. If $X\prec_{LMS}Y$ (see Section 2.3), then all the
suffixes of $T^{i}$ starting at positions in $\mathcal{X}$ are
lexicographically greater than the suffixes starting at positions in
$\mathcal{Y}$.

######  Proof 3.2.

Assume first that $X$ is not a prefix of $Y$ (and vice versa). We compare the
sequences of these strings from left to right until we find a mismatching
position $u$ (i.e., $X[u]\neq Y[u]$). We know that symbols $X[u]$ and $Y[u]$
define the lexicographical order of the suffixes in $\mathcal{X}$ relative to
the suffixes in $\mathcal{Y}$. In the other scenario, when one string is a
prefix of the other, we can not use this mechanism as we will not find a
mismatching position $X[u]\neq Y[u]$. For this case, we resort to the symbol
types of Section 2.3. We assume for this proof that $X$ is a prefix of $Y$,
but the other way is equivalent. We know that $X[x]$ and $Y[x]$ have different
types. $X[x]$ is LMS type because $X$ is a suffix of an LMS substring. On the
other hand, $Y[x]$ is L type because if it were S type, then it would also be
LMS type, and thus $Y[1,x]$ would be an occurrence for $X$. This observation
is due to $Y[x-1]=X[x-1]$ is L type. Given the types of $X[x]$ and $Y[x]$, the
occurrences of $X$ in $\mathcal{X}$ are always followed in $T^{i}$ by symbols
that are greater than $Y[x+1]$, meaning that the suffixes of $T^{i}$ starting
at positions in $\mathcal{X}$ are lexicographically greater than the suffixes
starting at positions in $\mathcal{Y}$. This observation does not hold when
$X$ or $Y$ have length one: $X[x]$ equals $Y[1]$ and both are LMS type, so
there is no enough information to decide the lexicographical order of the
suffixes in $\mathcal{X}$ and $\mathcal{Y}$.

The consequence of Lemma 3.1 is that the suffixes of length $>1$ in $D^{i}$
induce a partition over $S\\!A^{i}$ (the suffix array of $T^{i}$):

######  Lemma 3.3.

Let $\mathcal{S}=\\{S_{1},S_{2},\ldots,S_{k}\\}$ be the set of strings of
length $>1$ that occur as suffixes in the phrases of $D^{i}$. Additionally,
let $\mathcal{O}=\\{O_{1},O_{2},\ldots,O_{k}\\}$ be the set of occurrences in
$T^{i}$ for the strings in $\mathcal{S}$. For every $S_{u}\in\mathcal{S}$, its
associated list $O_{u}\in\mathcal{O}$ stores each position $j$ such that
$T^{i}[j,j+|S_{j}|-1]$ is an occurrence of $S_{u}$ and
$T^{i}[j-j^{\prime},j+|S_{j}|-1]$, with $j^{\prime}\geq 0$, is an LMS
substring. It holds that $\mathcal{O}$ induces a partition over the suffix
array of $T^{i}$ ($S\\!A^{i}$) as the lexicographical sorting places the
elements of each $O_{u}\in\mathcal{O}$ in a consecutive range of $S\\!A^{i}$.

######  Proof 3.4.

We demonstrate the lemma by showing that the lexicographical sorting does not
interleave suffixes of $T^{i}$ in $S\\!A^{i}$ that belong to different lists
of $\mathcal{O}$. Assume a string $S_{u}\in\mathcal{S}$, associated with the
list $O_{u}\in\mathcal{O}$, is a prefix in another string
$S_{u^{\prime}}\in\mathcal{S}$, which in turn is associated with the list
$O_{u^{\prime}}\in\mathcal{O}$. Even though we do not know the symbols that
occur to the right of $S_{u}$ in its occurrences of $O_{u}$, we do know that
both $S_{u}$ and $S_{u^{\prime}}$ are suffixes of LMS substrings, and by Lemma
3.1, we know that all the suffixes of $T^{i}$ in $O_{u}$ are lexicographically
greater than the suffixes in $O_{u^{\prime}}$. Hence, the interleaving of
suffixes in $S\\!A^{i}$ from different lists of $\mathcal{O}$ is not possible,
even if $\mathcal{S}$ is not a prefix-free set.

Lemma 3.3 gives us a simple way to construct the preliminary BWT for $T^{i}$
($pBWT^{i}$). We consider for the moment $pBWT^{i}$ to be a vector of lists to
simplify the explanations. We first sort the strings of $\mathcal{S}$ in
$\prec_{LMS}$ order. Then, for every $oth$ string $S\in\mathcal{S}$ in
$\prec_{LMS}$ order, we insert in the list $pBWT^{i}[o]$ the symbols that
occur to the left of $S$ in $D^{i}$. There are three cases to consider for
this task:

######  Lemma 3.5.

Let $S\in\mathcal{S}$ be the string with $\prec_{LMS}$ order $o$ among the
other strings in $\mathcal{S}$. If $S$ is left-maximal in $D^{i}$, then the
list $pBWT^{i}[o]$ contains more than one distinct symbol, and it is not
possible to decide the relative order of those symbols with the information of
$D^{i}$.

######  Proof 3.6.

Let $X$ and $Y$ be two phrases of $D^{i}$ where $S$ occurs as a suffix. Assume
the left symbol of $S$ in $X$ is $x\in\Sigma^{i}$ and the left symbol in $Y$
is $y\in\Sigma^{i}$. In this scenario, the relative order of $x$ and $y$ is
not decided by $S$, but for the sequences that occur to the right of $X$ and
$Y$ in $T^{i}$. However, those sequences are not accessible directly from
$D^{i}$. Hence, it is not possible to decide the order of $x$ and $y$ in
$pBWT^{i}[o]$.

######  Lemma 3.7.

Consider the string $S\in\mathcal{S}$ of Lemma 3.5. When $S$ occurs as a non-
proper suffix in a phrase $F\in D^{i}$, it is not possible to complete the
sequence of symbols for $pBWT^{i}[o]$.

######  Proof 3.8.

The symbols that occur to the left of $S$ in $T^{i}$ are stored in the LMS
substrings that precede $F$ in $T^{i}$. However, it is not possible to know
from $D^{i}$ which are those substrings.

We now describe the information of $pBWT^{i}$ that we can extract from
$D^{i}$:

######  Lemma 3.9.

Let $S\in\mathcal{S}$ be the string of Lemma 3.5. Additionally, let
$O\in\mathcal{O}$ be the list of occurrences of $S$ in $T^{i}$ as described in
Lemma 3.3. If all the suffixes of $T^{i}$ in $O$ are preceded by the same
symbol $s\in\Sigma^{i}$ (i.e., $S$ is not left-maximal), then
$pBWT^{i}[o]=(s,l)$ is an equal-symbol run of length $l=|O|$, where $o$ is the
$\prec_{LMS}$ order of $S$ in $\mathcal{S}$.

######  Proof 3.10.

By Lemma 3.3, we know that the suffixes of $T^{i}$ in $O$ are prefixed by $S$,
and that they form a consecutive range $S\\!A^{i}[j,j^{\prime}]$.
Additionally, the symbols that occur to the left of the suffixes in
$S\\!A^{i}[j,j^{\prime}]$ are those for the list of $pBWT^{i}[o]$. However, we
still have not resolved the relative order of the suffixes in
$S\\!A^{i}[j,j^{\prime}]$, so (in theory) we do not know how rearrange the
symbols in $pBWT^{i}[o]$. The suffixes of $T^{i}$ in $O$ are preceded by the
same symbol $s$, so it is no necessary to further sort
$S\\!A^{i}[j,j^{\prime}]$ because the outcome for $pBWT^{i}[o]$ will be always
an equal-symbol run for $s$ of length $l=|O|$.

The problem is that we do not store $O$, so we do not know value for $l$ in
$(s,l)$. Nevertheless, we do have the frequencies of the phrases in $D^{i}$,
in the vector $N^{i}$. In this way, we can compute $l$ by summing the
frequencies in $N^{i}$ for the phrases of $D^{i}$ where $S$ occurs as a
suffix.

Now that we have covered all the theoretical aspects of the parsing phrase, we
proceed to describe our procedure to build $pBWT^{i}$.

####  3.3.1 Constructing the Preliminary BWT for the Parsing Round

The computation of $pBWT^{i}$ starts with the construction of a _generalized_
suffix array $S\\!A_{D^{i}}$ for $D^{i}$. We say $S\\!A_{D^{i}}$ is
generalized because it only considers the suffixes of the dictionary phrases.
If a string $S\in\mathcal{S}$ appears as a suffix in two or more phrases,
those occurrences maintain in $S\\!A_{D^{i}}$ the relative order in which
their enclosing phrases appear in $D^{i}$. In practice, the values we store in
$SA_{D^{i}}$ are the positions in $R^{i}$, the vector storing the concatenated
phrases of $D^{i}$ (see the encoding of $D^{i}$ in Subsection 3.3).

We compute $S\\!A_{D^{i}}$ using a modified version of the ISS method
mentioned in Subsection 2.3. The first difference is that, in step one, we
insert in $S\\!A_{D^{i}}$ the position in $R^{i}$ of the last symbol of each
phrase. Put it another way, suppose $R^{i}[j]$, with $L^{i}[j+1]=1$, is the
last symbol of a phrase $F$, then we insert $j$ in the right-most available
cell in the bucket $R^{i}[j]$ of $S\\!A_{D^{i}}$. The step one in the original
ISS puts LMS-type symbols at the end of the buckets. In our case, the last
symbol of a phrase is, by definition, LMS type in $T^{i}$, so the operation is
homologous. The second difference of our ISS variation is that, during step
two and three, we skip each position $S\\!A_{D^{i}}[u]$ representing the start
of a phrase ($L^{i}[S\\!A_{D^{i}}[u]]=1$) as they do not induce suffixes.

The next step is to scan $S\\!A_{D^{i}}$ from left to right to compute
$pBWT^{i}$. From now on, we consider $pBWT^{i}$ to be a run-length compressed
vector instead of a vector of lists. As we move throughout the suffix array,
we search for every range $S\\!A_{D^{i}}[j,j^{\prime}]$, with
$j^{\prime}-j+1\geq 1$, that encode suffixes with the same sequence222In
practice, we compute every distinct range $S\\!A_{D^{i}}[j,j^{\prime}]$ during
the construction of the suffix array. We reserve the least significant bit in
the cells of $S\\!A_{D^{i}}$ to mark every position $S\\!A_{D^{i}}[j]$. We
flag these positions during the execution of our modified version of ISS
(Section 2.3).. Nevertheless, we consider only the ranges that either
represent suffixes of length $>1$ or suffixes of length $1$ that expand to
suffixes of $\mathcal{T}$. Recall that the left-most symbol of an LMS
substring is the same as the right-most symbol of the LMS substring that
precedes it. Hence, considering all the suffixes in $D^{i}$ will produce a
redundant (and incorrect) BWT. The only exception to this rule are the LMS
substrings at the beginning of the strings of $\mathcal{T}$ as they do not
share a symbol with the LMS substring to their left. This kind of substrings
only appear when we cross from $T_{u+1}$ to $T_{u}$ in $T^{i}$, with
$T_{u},T_{u+1}\in\mathcal{T}$. We can detect this situation using $B^{i}$, the
bit vector marking the symbols in $\Sigma^{i}$ that expand to suffixes of
strings in $\mathcal{T}$ (see Section 3.3).

We define the length of $S\\!A_{D^{i}}[j,j^{\prime}]$ as
$l=\sum_{u=j}^{j^{\prime}}N^{i}[\mathsf{rank}(L^{i},S\\!A_{D^{i}}[u])]$. This
value is the sum of the frequencies of the phrases where the suffixes in
$S\\!A_{D^{i}}[j,j^{\prime}]$ occur.

If all the suffixes of $S\\!A_{D^{i}}[j,j^{\prime}]$ are followed by the same
symbol $s\in\Sigma^{i}$, we append $(s,l)$ to $pBWT^{i}$ (see Lemma 3.9).
Otherwise we append $(\texttt{*},l)$. The symbol * represents an empty entry
and it is out of $\Sigma^{i}$. We will resolve $(\texttt{*},l)$ in the next
phase of $\mathsf{grlBWT}$ (see Lemmas 3.5 and 3.7). After scanning
$S\\!A_{D^{i}}$, we store $pBWT^{i}$ into disk and discard $S\\!A_{D^{i}}$.

####  3.3.2 Grammar Compression and Next Parsing Round

Once we finish constructing $pBWT^{i}$, the next step in the parsing round $i$
is to store $D^{i}$ in a compact form to use it later during the induction
phase of $\mathsf{grlBWT}$. We first explain why we need $D^{i}$ during the
induction phase and then describe the format we choose to encode it.

Broadly speaking, the induction process consists of scanning $BWT^{i+1}$ from
left to right, mapping every symbol $BWT^{i+1}[j]\in\Sigma^{i+1}$ back to the
phrase $F\in D^{i}$ from which it originated, and then checking which of the
proper suffixes of $F$ produced empty entries in $pBWT^{i}$ (see Lemmas 3.5
and 3.7). Assume the suffix $F[u..]=S\in\mathcal{S}$ produced an empty entry,
then we append $F[u-1]$ in the BWT range associated with $S$ (see Lemma 3.3).

The process described above requires $D^{i}$ and a mechanism to map the left-
maximal suffixes in $D^{i}$ back to the empty entries they produce in
$pBWT^{i}$. We solve the problem by encoding $D^{i}$ with a representation
that is similar to grammar compression (Section 2.2).

We start by discarding $N^{i}$ and the $\textsf{rank}_{1}$ data structure, as
they are no longer necessary (see the current encoding of $D^{i}$ in Section
3.3). For our method to work, we also need each string $S\in\mathcal{S}$
associated with an empty entry $(\texttt{*},l)$ of $pBWT^{i}$ to be a member
of $D^{i}$. This property might not hold when the string $S$ that produced an
empty entry meets Lemma 3.5. The problem arises if $S$ always appears as a
proper suffix in $D^{i}$, not as a full phrase. If that is the case, we
_create_ 333It means we append the sequence of $F$ at the end of $R^{i}$ and
expand $L^{i}$ accordingly. a new independent entry for $S$ in $D^{i}$.

After expanding $D^{i}$, we create a hash table $M^{i}$ in which we insert
each phrase $F\in D^{i}$ occurring as a left-maximal suffix. If $F$ has
$\prec_{LMS}$ rank $b$ in $D^{i}$, then we insert the pair $(F,b)$ into
$M^{i}$, where $F$ is the key and $b$ is the value. Once we construct $M^{i}$,
we reorder the way in which the phrases of $D^{i}$ are concatenated in $R^{i}$
according to their $\prec_{LMS}$ ranks.

The next step consists of _compressing_ $D^{i}$. We scan $R^{i}$ from left to
right, and for every $F=R^{i}[j,j^{\prime}]\in D^{i}$, with $L^{i}[j]=1$ and
$L^{i}[j^{\prime}+1]=1$, we search for the longest proper suffix $F[u..]$ that
exists in $M^{i}$ as a key. If such key exists, then we replace $F$ with
$R[j,j+1]=F[u-1]{\cdot}b^{\prime}$, where $b^{\prime}$ is the value associated
with $F[u..]$ in $M^{i}$. If no proper suffix of $F$ exists in $M^{i}$, then
we replace $F$ as $R[j,j+1]=\texttt{*}{\cdot}F[|F|-1]$, where * is a dummy
symbol. After updating the sequence of $F$, we mark the symbols in
$R[j+2,j^{\prime}]$ as discarded if $j^{\prime}-j+1=|F|\geq 2$. When we finish
the scan of the dictionary, we left-contract $R^{i}$ by removing the discarded
symbols. This process reduces the phrases in $D^{i}$ to strings of length two,
so the vector $L^{i}$ is no longer necessary.

Now we explain the rationale of our encoding. We develop our argument as a
chain of implications. Consider again the phrase $F$, which we replaced with
the sequence $F[u-1]{\cdot}b^{\prime}$. We obtained
$b^{\prime}\in\Sigma^{i+1}$ when we performed a lookup operation of $S=F[u..]$
in $M^{i}$ during the compression of $D^{i}$. The value $b^{\prime}$ that the
lookup returned is the $\prec_{LMS}$ order of $S$ in $D^{i}$. The membership
of $S$ to the keys of $M^{i}$ implies that $S$ appears as a left-maximal
suffix in $D^{i}$, which in turn implies that $S$ is a full phrase in $D^{i}$
too (we enforced this property when we expanded the dictionary). Additionally,
the left-maximal condition of $S$ implies that there were at least two suffix
occurrences of $S$ preceded by different symbols. This is why $S$ produces an
empty entry in $pBWT^{i}$. Now, recall that we sorted the phrases of $D^{i}$
in $\prec_{LMS}$ order in $R^{i}$. Therefore, if we want to access $S$, we
have to go to the substring $R^{i}[2b^{\prime}-1,2b^{\prime}]$. This substring
does not encode the full sequence of $S$, but its longest left-maximal suffix
$R^{i}[2b]\in\Sigma^{i+1}$ (which is also a left-maximal suffix of $F$) along
with the left-context symbol for that suffix ($R^{i}[2b-1]\in\Sigma^{i}$).
Recursively, the longest left-maximal symbol of $S$ is not a sequence either,
but a pointer to another position of $R^{i}$. We access this nested left-
maximal suffix by setting $b^{\prime}=R^{i}[2b^{\prime}]$ and updating the
values $R^{i}[2b^{\prime}-1,2b^{\prime}]$. We continue applying this idea
until we reach a range $R^{i}[2b^{\prime}-1,2b^{\prime}]$ where
$R^{i}[2b^{\prime}]\in\Sigma^{i}$, which implies that we reached the last
suffix of $F$. This last range will store the sequence
$\texttt{*}{\cdot}F[|F|-1]$. Notice that the right symbol in this case is not
the last symbol of $F$ but its left context $F[|F|-1]$. This is because the
LMS substrings overlap by one character in $T^{i}$, so $F[|F|]$ is redundant
as it also appears as a prefix in another phrase. However, we need $F[|F|-1]$,
as we will append it to one of the empty entries of $pBWT^{i}$ in the next
phase of $\mathsf{grlBWT}$. The only exception to this rule is when $F$
expands to a suffix of a string in $\mathcal{T}$. In that case, we store
$F[|F|]$ instead of $F[|F|-1]$ in $R[2b^{\prime}-1,2b^{\prime}]$ as
$F[|F|]=R^{i}[2b^{\prime}]$ is not a prefix in any other phrase. On the other
hand, we need the dummy symbol * to maintain the invariant that all the
phrases encoded in $R^{i}$ have length two. Once we finish the compression, we
store $R^{i}$ on disk. From now on, we use $D^{i}$ to refer to $R^{i}$.

Figure 1: Parsing round $i=1$ of $\mathsf{grlBWT}$ for the collection
$\mathcal{T}=\\{\texttt{gtacc},\texttt{gtaatagtacc}\\}$, with
$T=\texttt{gtacc\$gtaatagtacc\$}$. (A) Construction of the dictionary
$D^{1}=\\{R^{1},L^{1},N^{1}\\}$ from $T^{1}=T$. The sequence in gray below
$T^{1}$ stores the symbol types: L-type is $L$, S-type is $S$ and LMS-type is
$S^{*}$. The dashed vertical lines in $R^{1}$ mark the boundaries between
dictionary phrases. (B) Constructing $pBWT^{1}$ from the suffix array
$S\\!A_{D^{i}}$ of $D^{i}$. The ranges of $SA_{D^{1}}$ enclosed by dashed
boxes produce empty entries in $pBWT^{i}$. Notice we do not use the entries
$SA_{D^{i}}[5,7]=3,7,11$ for computing $pBWT^{1}$ as their corresponding
symbols are redundant. For instance, consider $R^{1}[3]=\texttt{a}$, which is
a suffix in $R^{1}[1,3]=\texttt{gta}$. The last symbol in the occurrences of
gta in $T^{1}$ overlaps the first symbol of acc$ or aata. Therefore, the
symbol $R^{1}[3]=\texttt{a}$ of gta is covered by acc$ or aata in $pBWT^{1}$.
We represent $pBWT^{1}$ as a sequence of equal-symbol runs. The upper row
depicts the run symbols while the lower row show the run lengths. (C)
Compressing $D^{1}$. The strings at the top are the dictionary phrases sorted
in $\prec_{LMS}$ order. We add the string ta to the dictionary as it appears
as a left-maximal suffix in $D^{i}$, and hence, produces an empty entry in
$pBWT^{1}$. The sequences in gray are the longest left-maximal suffixes of the
phrases. The underlined symbols are the left contexts of those suffixes.
Notice we compress $F=\texttt{acc\$}$ directly to *$ as $F$ does not have a
left-maximal suffix. Besides, as $F$ does not overlap other phrases in $T^{1}$
(because of $), we store $F[|F|]=\texttt{\$}$ instead of
$F[|F|-1]=\texttt{c}$. A different situation occurs with $S=\texttt{ta}$. This
phrase does not have left-maximal suffixes of length $>1$. However, in this
case, we store $S[|S|-1]=\texttt{t}$ instead of $S[|S|]=\texttt{a}$ as $S$
overlaps other phrases in $T^{1}$. (D) The parse $T^{2}$ we obtain by
replacing the phrases of $T^{1}$ with their $\prec_{LMS}$ orders in $D^{i}$.
The gray symbols expand to suffixes of strings in $\mathcal{T}$.

The final step for the parsing round $i$ is to create the new text $T^{i+1}$.
We first reload from disk the hash table $H^{i}$ we produced at the beginning
of the iteration, and replace its values with the keys’ $\prec_{LMS}$ orders.
More specifically, if a key $F\in D^{i}$ of $H^{i}$ has $\prec_{LMS}$ order
$b$ among the other strings that produced empty entries in $pBWT^{i}$, then we
update the value of $F$ in $H^{i}$ to $b$. Notice that the strings we stored
as keys in $H^{i}$ are not the same as those we have now in $D^{i}$ because we
compress them. Therefore, we can not lookup the phrases of $D^{i}$ in the keys
of $H^{i}$ to update the hash table values. Still, we can overcome this
problem if we modify $H^{i}$ after sorting $D^{i}$ in $\prec_{LMS}$ order but
before we compress it.

Once we update $H^{i}$, we construct $T^{i+1}$ by scanning $T^{i}$ again and
replacing the LMS substrings with their associated values in $H^{i}$. If
$T^{i+1}$ has length $k$ (the number of strings in $\mathcal{T}$), then we
stop the parsing phase as all the strings in $\mathcal{T}$ are now compressed
to one symbol. An example of the parsing step is depicted in Figure 1.

###  3.4 The Induction Phase

The induction phase starts with the computation of $BWT^{h}$, the BCR BWT for
the text $T^{h}$ of the last parsing round $h$. This step is trivial as each
symbol in $T^{h}$ encodes a full string of $\mathcal{T}$ (see the ending
condition of the parsing phase). Hence, the left context of every symbol is
the symbol itself. BCR BWT maintains the relative order of the strings in
$\mathcal{T}$ (see Section 2.1), so $BWT^{h}$ is $T^{h}$ itself.

We now describe the steps we perform during every induction step $i<h$. In
this case, assume we receive as input (i) $BWT^{i+1}$ from the previous phase,
(ii) $D^{i}$, and (iii) $pBWT^{i}$. Before explaining our procedure, we
describe some important properties of $BWT^{i+1}$.

######  Lemma 3.11.

Let $BWT^{i+1}[j]$ and $BWT^{i+1}[j^{\prime}]$ be two symbols at different
positions $j$ and $j^{\prime}$, with $j<j^{\prime}$, whose mapping phrases in
$D^{i}$ are $F$ and $F^{\prime}$, respectively. Also, let the proper suffixes
$F[u..]=F^{\prime}[u^{\prime}..]=S\in\mathcal{S}$ (see Lemma 3.3 for the
description of $\mathcal{S}$) be two occurrences in $T^{i}$ of a string $S$
that appears as a left-maximal suffix in $D^{i}$. The suffix of $T^{i}$
prefixed by $F[u..]$ precedes in $S\\!A^{i}$ (the suffix array of $T^{i}$) the
suffix of $T^{i}$ prefixed by $F[u^{\prime}..]$.

######  Proof 3.12.

As $F[u..]$ and $F^{\prime}[u^{\prime}..]$ are equal, their relative orders
are decided by the right contexts in $T^{i}$ of the occurrences $BWT^{i+1}[j]$
and $BWT^{i+1}[j^{\prime}]$ of $F$ and $F^{\prime}$ (respectively). By
induction, we know that $BWT^{i+1}$ is complete, and as $BWT^{i+1}[j]$
precedes $BWT^{i+1}[j^{\prime}]$ in the BWT, the right context of $F[u..]$ has
a smaller order in $S\\!A^{i}$ than the right context of
$F^{\prime}[u^{\prime}..]$.

If we generalize Lemma 3.11 to $x\geq 1$ occurrences of $S$, then we can use
the following lemma to compute the sequence for the empty entry of $pBWT^{i}$
generated by $S$:

######  Lemma 3.13.

Let $S$ be a string of $\mathcal{S}$. Additionally, let
$J=\\{j_{1},j_{2},\ldots,j_{x}\\}$ be a list of strictly increasing positions
of $BWT^{i+1}$. Every $BWT^{i+1}[j_{o}]$, with $j_{o}\in J$, is a symbol
$b\in\Sigma^{i+1}$ generated from a phrase $F\in D^{i}$ where $S=F[u..]$
occurs as a proper suffix. The symbols of $B^{i+1}$ referenced by $J$ are not
necessarily equal, and hence, their associated phrases in $D^{i}$ are not
necessarily the same. However, these phrases of $D^{i}$ are all suffixed by
$S$. Assume we scan $J$ from left to right, and for every $j_{o}$, we extract
the symbol $F[u-1]\in\Sigma^{i}$ that precedes $S$ and append it to a list
$L_{S}$. The resulting list $L_{S}\in\Sigma^{i*}$ has the same sequence of
symbols as the $BWT^{i}$ range that maps the block for $S$ in the partition of
$\mathcal{S}$.

######  Proof 3.14.

Because of Lemma 3.11, the suffix of $T^{i}$ prefixed by the occurrence
$BWT^{i+1}[j_{o}]$ of $S$ precedes the suffix of $T^{i}$ prefixed by the
occurrence $BWT^{i+1}[j_{o+1}]$. This property holds for every $j_{o}$, with
$o\in[1,x-1]$. Hence, the suffixes of $T^{i}$ prefixed by $S$ are already
sorted $J$.

Our compressed representation for $D^{i}$ (see Section 3.3.2) has precisely
the information we need to construct $L_{S}$ as described in the procedure of
Lemma 3.13. Still, the idea only works when $S$ _always_ appears as a proper
suffix in the phrases of $D^{i}$. When $S$ matches a full phrase (see Lemma
3.7), there is no left-context symbol for $S$ we can extract from $D^{i}$.
Nevertheless, there is only one phrase $F\in D^{i}$ where $S$ can be a non-
proper suffix, and because $S$ comes from $BWT^{i+1}$, $F$ has to be an LMS
substring in $T^{i}$. This observation implies that $F$ maps to a symbol
$b\in\Sigma^{i+1}$ in $T^{i+1}$. Hence, we can extract the left-context
symbols of $S$ from the range in $BWT^{i+1}$ that corresponds to the $bth$
bucket of $S\\!A^{i+1}$. We explain how to carry out this process in the next
subsection.

###  3.5 The Induction Algorithm

Let $p$ be the sum of the lengths in the empty entries of $pBWT^{i}$. These
lengths correspond to the second field in the run-length representation of
$pBWT^{i}$. We start the induction by creating a vector $P^{i}[1,p]$, which we
logically divide into $\sigma^{i+1}$ buckets (recall that $\sigma^{i+1}$
matches the number of empty entries in $pBWT^{i}$). Every bucket $b$ of
$P^{i}$ will be of size $l_{b}$, the length of the $bth$ empty entry (from
left to right) of $pBWT^{i}$. Subsequently, we perform a scan over $BWT^{i+1}$
from left to right. For each symbol $BWT^{i+1}[j]=b\in\sigma^{i+1}$, we first
check if its associated LMS substring $F\in D^{i}$ (the string from which we
obtain the symbol $b$ during the parsing round $i$) exists as a suffix in
other phrases of $D^{i}$. This information is already encoded in a bit vector
$V^{i}[1,\sigma^{i+1}]$ we constructed during the parsing round $i$. When $F$
occurs as an LMS substring and as a proper suffix in other dictionary phrases
($V^{i}[b]=1$), we append a dummy symbol in the bucket $b$ of $P$. This is the
situation we described at the end of the previous subsection. After processing
$b$, we decompress the left-maximal suffixes of its phrase $F$ from the
compressed representation of $D^{i}$.

The decompression of $F$ begins by accessing the range $D^{i}[2b-1,2b]$ (see
Section 3.3.2). If $o=D^{i}[2b]$ belongs to $\Sigma^{i+1}$, then the symbol
$o$ encodes a string $F[u..]=S$, with $u>1$, whose sequence is a left-maximal
suffix in $D^{i}$. During the parsing phase of $\mathsf{grlBWT}$, we inserted
$S$ to $D^{i}$ as an independent string as it yields an empty entry for
$pBWT^{i}$ (see Lemma 3.5). The order of $S$ in $D^{i}$ is precisely $o$, its
$\prec_{LMS}$ rank among the other phrases of $D^{i}$. On the other hand, the
left-context symbol of $S$ is $D^{i}[2b-1]\in\Sigma^{i}$. With this
information, we apply Lemma 3.13 by appending the symbol $D^{i}[2b-1]$ to the
bucket $o$ of $P$. Then, we move to the next left-maximal suffix of $F$ by
setting $b=o$ and updating the range $D^{i}[2b-1,2b]$.

The decompression of $F$ stops when $D^{i}[2s]$ belongs to $\Sigma^{i}$, which
means we reach the last symbol of $F$. For the moment, we do not know for
which phrase of $D^{i}$ $D^{i}[2s]$ is its left context. Hence, we set
$BWT^{i+1}[j]=D^{i}[2s]$ and leave this position on hold to process it later.
After finishing the scan of $BWT^{i+1}$, its symbols are now over the alphabet
$\Sigma^{i}$. These values are the ones we have to insert in the dummy
positions of $P^{i}$. Notice that the entries in $P^{i}$ and $BWT^{i+1}$ are
already sorted by their right contexts in $T^{i}$. Hence, the completion of
the dummy positions reduces to a merge of two sorted lists.

The last step in the induction round $i$ consists of merging $pBWT^{i}$,
$BWT^{i+1}$ and $P^{i}$ in $BWT^{i}$. We scan $pBWT^{i}$ and we append its
entries to $BWT^{i}$ as long as they are not empty. Then, when we reach an
empty entry $(\texttt{*},l)$, we proceed as follows: assume the current pair
$(\texttt{*},l)$ is the $bth$ empty entry of $pBWT^{i}$. Then, we check if the
phrase $F\in D^{i}$ that produced this entry (the one with $\prec_{LMS}$ order
$b$) only occurs as a full LMS substring in $T^{i}$ ($V^{i}[b]=0$). If that is
the case, we append the next $l$ symbols of $BWT^{i+1}$ into $BWT^{i}$. On the
other hand, when $F$ occurs as an LMS substring, but also as a proper suffix
in other phrases of $D^{i}$ ($V^{i}[b]=1$), the next $l$ symbols of $BWT^{i}$
are a mix of entries from the bucket $b$ of $P^{i}$ and $BWT^{i+1}$. We append
symbols from the bucket $b$ of $P^{i}$ as long as they are not dummy. When we
reach a dummy symbol in $P^{i}$, we change the list, and append the next $x$
symbols of $BWT^{i+1}$ into $BWT^{i}$, where $x$ is the number of consecutive
dummy symbols we saw in $P^{i}$. Once we process the $x$ entries of
$BWT^{i+1}$, we go back to $P$ and continue back and forth between $P$ and
$BWT^{i}$ until we process all the symbols in the bucket $b$ of $P$.

The last case we have to cover for the merge is when $F$ always occurs as a
proper suffix in $D^{i}$ (i.e., it is not an LMS substring of $T^{i}$). This
situation is simple as we marked $F$ in $V^{i}$ ($V^{i}[b]=1$). Hence, we just
copy the content of the bucket $b$ of $P$ into $BWT^{i}$. Notice this bucket
will not have dummy entries as $b$ does not appear in $B^{i+1}$ as a symbol.
We obtain occurrences of $b$ while decompressing other phrases of $D^{i}$
whose $\prec_{LMS}$ ranks do appear as symbols in $B^{i+1}$, and where $F$ is
a proper left-maximal suffix. Once we complete all the induction rounds, the
final BCR BWT is in $BWT^{i}$.

####  3.5.1 Speeding up the Induction with Compression

If we run-length encode $BWT^{i+1}$, the induction becomes more efficient.
Every position $BWT^{i+1}[j]$ is not a symbol $b$, but a pair $(l,b)$ that
represents $l$ consecutive copies of $b$. Thus, instead of decompressing $l$
times the left-maximal suffixes of the phrase associated with $b$, we
decompress them only once, and copy the result $l$ times to the different
buckets of $P^{i}$.

Maintaining $P^{i}$ as a run-length encoded sequence also improves efficiency.
A compact representation of $P^{i}$ reduces the working memory, which in turn
reduces the number of cache misses. The only problem with this idea is that we
do not know before the induction how many equal-symbols runs $P^{i}$ will
have. There are two solutions to this problem. First, we could represent
$P^{i}$ as a dynamic vector [2]. Its initial size would be $\sigma^{i+1}$, and
then we expand the buckets as we insert new equal-symbol runs into them. The
second option is to perform a preliminary scan of $BWT^{i+1}$ to compute the
size of the run-length compressed version of $P^{i}$, then we scan $BWT^{i+1}$
again to perform the induction.

###  3.6 Complexity of our Method

We show that the construction of the BWT remains linear, even though we
perform compression during the intermediate steps.

######  Theorem 3.15.

Let $\mathcal{T}=\\{T_{1},T_{2},\ldots,T_{k}\\}$ be a collection with $k$
strings and $n$ symbols. The algorithm $\mathsf{grlBWT}$ constructs the BCR
BWT of $\mathcal{T}$ in expected $O(n)$ time and requires $O(n)$ bits of
working space.

######  Proof 3.16.

The complexities in the theorem were already proved for the linear-time
algorithms that construct the suffix array [30] and the BWT [34] using ISS. We
show that these complexities are not altered by our compression scheme. Let
$n^{i}=|T^{i}|$ be the length of the input text we receive at parsing round
$i$. Hashing the dictionary phrases from $T^{i}$ runs in $O(n^{i})$ expected
time and requires $O(n^{i})$ bits of space. The construction of
$S\\!A_{D^{i}}$ runs in $O(n^{i})$ time and space as we use ISS to build it,
and the number of symbols in $D^{i}$ is never greater than $n^{i}$. The extra
steps of the parsing round only require a constant number of linear scans over
$S\\!A_{D^{i}}$. During the induction phase, we only perform linear scans over
$BWT^{i}$ and $pBWT^{i}$. We still have the cost of accessing the left-maximal
suffixes of $D^{i}$ when we scan $BWT^{i+1}$ during the induction phase.
However, our simple compressed representation for $D^{i}$ (Section 3.3.2)
supports random access in $O(1)$ time to the symbols, and the number of left-
maximal suffixes we visit during the scan of $BWT^{i+1}$ is no more than
$n^{i}$. In every parsing round, the size of $T^{i+1}$ is at most half the
size of $T^{i}$, so the sum of the text lengths $n^{1},n^{2},\ldots,n^{h}$ is
$O(n)$ (see Nong et al. [30]). This property also implies that grlBWT never
visits more than $O(n)$ left-maximal suffixes during its induction phase.

Dataset | $\sigma$ | Number of strings | Max. length | Avg. length | $n$ | $n/r$
---|---|---|---|---|---|---
ILL1 | 5 | 84,006,956 | 151 | 151 | 12,769,057,312 | 3.18
ILL2 | 5 | 160,285,798 | 151 | 151 | 24,363,441,296 | 4.07
ILL3 | 5 | 235,805,550 | 151 | 151 | 35,842,443,600 | 4.67
ILL4 | 5 | 305,931,740 | 151 | 151 | 46,501,624,480 | 5.03
ILL5 | 5 | 377,453,488 | 151 | 151 | 57,372,930,176 | 5.33
HGA05 | 16 | 334,065 | 248,956,422 | 42,715 | 14,269,998,434 | 4.82
HGA10 | 16 | 759,341 | 250,522,664 | 39,025 | 29,634,170,092 | 8.76
HGA15 | 16 | 835,485 | 250,522,664 | 53,918 | 45,048,695,199 | 12.02
HGA20 | 16 | 874,235 | 250,522,664 | 68,650 | 60,017,146,889 | 15.67
HGA25 | 16 | 899,424 | 250,522,664 | 83,447 | 75,055,723,570 | 19.42
Table 1: Datasets. The upper rows are the Illumina reads while the lower rows
are the human genomes. Columns four and five are the minimum and average
string length (respectively) in the collection. The value for $r$ is the
number of equal-symbol runs in the BCR BWT of the collection.

##  4 Experiments

We implemented $\mathsf{grlBWT}$ as a C++ tool, also called grlBWT. This
software uses the SDSL-lite library [11] to operate with bit vectors and rank
data structures. Our source code is available at
https://github.com/ddiazdom/grlBWT. We compared the performance of grlBWT
against other tools that compute BWTs for string collections:

  * •

ropebwt2444https://github.com/lh3/ropebwt2 : a variation of the original BCR
algorithm of Bauer et al. [1] that uses rope data structures [2]. This method
is described in Heng Lee [20].

  * •

pfp-eBWT555https://github.com/davidecenzato/PFP-eBWT : the eBWT algorithm of
Boucher et al. [4] that builds on PFP and ISS.

  * •

BCR_LCP_GSA666https://github.com/giovannarosone/BCR_LCP_GSA : the current
implementation of the semi-external BCR algorithm [1].

  * •

egap777https://github.com/felipelouza/egap: a semi-external algorithm of Edigi
et al. [8] that builds the BCR BWT.

  * •

gsufsort888https://github.com/felipelouza/gsufsort: an in-memory method
proposed by Louza et al. [23] that computes the BCR BWT and (optionally) other
data structures.

We also considered the tool bwt-lcp-em [3] for the experiments. Still, by
default it builds both the BWT and the LCP array, and there is no option to
turn off the LCP array, so we discarded it. We compiled all the tools
according to their authors’ description. For grlBWT, we used the compiler
flags -O3 -msse4.2 -funroll-loops.

We considered two common types of genomic data for the experiments: short
reads and assembled genomes. We downloaded five collections of Illumina reads
produced from different human genomes
999https://www.internationalgenome.org/data-portal/data-collection/hgdp. We
concatenated the strings so that our dataset 1 had one read collection,
dataset 2 had two collections, and so on. We named the files ILLX, where X is
the number of read collections concatenated. We also downloaded from
NCBI101010https://www.ncbi.nlm.nih.gov/assembly 25 collections of fully-
assembled human genomes. Like with the reads, we created the inputs for our
experiments so that every dataset has five more genomes than the previous one.
This setup aims to increase the repetitiveness as the collection size
increases. We named each file using the prefix HGA concatenated with the
number of genomes it had. The only prepossessing step we performed on the
genomes was to put every chromosome in one line and set all the characters to
upper case. All our inputs are described in Table 1.

We also investigated the effect of page cache [24, Ch. 16] in grlBWT. In every
parsing round $i$, we keep $T^{i}$ on disk, and linearly scan its file by
loading from disk to RAM one data chunk of $8$ MB at the time. Similarly, we
keep a buffer of $8$ MB in RAM for $T^{i+1}$, which we dump into disk every
time it gets full. We manipulate $BWT^{i}$ (respectively, $BWT^{i+1}$) in the
same way. We used the function posix_fadvise to turn off the page cache for
$T^{i}$, $T^{i+1}$, $BWT^{i}$, and $BWT^{i+1}$. Then we assessed the
performance of grlBWT on the assembled genomes using posix_advice and not
using it. We did not evaluate the effect of the page cache in the other tools.

We limited the RAM usage of egap to three times the input size. For
BCR_LCP_GSA, we turned off the construction of the data structures other than
the BCR BWT and left the memory parameters by default. In the case of
gsufsort, we used the flag \--bwt to build only the BWT. For ropebwt2, we set
the flag -L to indicate that the data was in one-sequence-per-line format, and
the flag -R to avoid considering the DNA reverse strands in the BWT. We ran
the experiments on the Illumina reads using one thread in all programs as not
all support multi-threading. For this purpose, we set the extra flag -P to
ropebwt2 to indicate single-thread execution. We tested the human genomes only
on ropebwt2, grlBWT and pfp-ebwt. By default, ropebwt2 uses four working
threads, so we set the same number of threads for grlBWT and pfp-ebwt. We did
not report results for pfp-ebwt with dataset ILL25 as the execution crashed.
We carried out the experiments on a machine with Debian 4.9, 736 GB of RAM,
and processor Intel(R) Xeon(R) Silver @ 2.10GHz, with 32 cores.

Figure 2: Memory peak usage (GBs) and elapsed time (in hours) for the Illumina
reads.

##  5 Results and Discussion

We summarize our experiments in Figures 2 and 3. The results we report for
grlBWT do not consider the use of posix_advice to turn off the page cache. In
Illumina reads, the fastest method was ropeBWT2, with a mean elapsed time of
4.14 hours. It is then followed by BCR_LCP_GSA, gsufsort, grlBWT, pfp-bwt, and
egap, with mean elapsed times of 9.58, 9.43, 10.05, 13.08, and 27.30 hours,
respectively (Figure 2B). We notice that grlBWT is competitive with
BCR_LCP_GSA and gsufsort. However, it gets slightly faster than them from
input ILL4 onward. We expected this behaviour since the largest datasets are
more repetitive.

Regarding the working space, the most efficient was BCR_LCP_GSA, with an
average memory peak of 5.73 GB. It is then followed by ropebwt2, with an
average memory peak of 26.64 GB. In both cases, the memory consumption
increases slowly with the input size. In the case of grlBWT, the memory peak
is more considerable; 42.20 GB on average, with a memory consumption that
grows faster than the previous methods (see Figure 2A). However, egap,
gsufsort, and pfp-ebwt are far more expensive, and their memory consumption
grow even faster. The tool egap uses 110.94 GBs on average. On the other hand,
pfp-ebwt and gsufsort have similar average memory peaks: 331.98 and 372.68
GBs, respectively.

Figure 3: Memory peak usage (GBs) and elapsed time (in hours) for the
assembled genomes.

In the repetitive datasets (human genomes), the results changed drastically
(see Figure 3). Our tool grlBWT outperformed ropebwt2 and pfp-ebwt in elapsed
time, with an average of 4.89 hours versus 20.95 and 9.55 hours of ropebwt2
and pfp-ebwt, respectively. As expected, the time for grlBWT grows smoothly
with the input size, while the time for ropeBWT grows fast. The time function
for pfp-ebwt also grows smoothly, but the results are still slower than those
of grlBWT (see Figure 3B). Regarding memory peak, ropebwt2 is the most
efficient tool, with a mean of 18.05 GB. Still, grlBWT obtained competitive
results, with an average peak of 20.38 GB. In this case, the memory
consumption growth in grlBWT is slightly steeper than in ropebwt2, but it
remains smooth. In contrast, pfp-ebwt has a more dramatic growth in memory
consumption, with an average memory peak of 156.74 GB (Figure 3A).

We observe that ropeBWT2 and pfp-ebwt performed well in one measure, but not
in both. In contrast, grlBWT maintained a low footprint for both measures,
elapsed time and memory consumption. This result demonstrates that our
strategy of keeping the intermediate data of the BWT algorithm in compressed
format works well when the text is repetitive.

Our experiments on the page cache showed there is an average slowdown of
$19\%$ in grlBWT when the cache is disabled with the function posix_advice.
This slowdown factor increases with the input size, being the lowest with
HGA05 ($12\%$) and the highest with HGA20 ($26\%$). This result is expected as
we are only using a static buffer of $8$ MB. A simple solution would be to set
a dynamic buffer that uses, say, $0.5\%$ of the input instead of the fixed $8$
MB.

##  6 Concluding Remarks

We introduced a method for building the BCR BWT that maintains the data of
intermediate stages in compressed form. The representation we chose not just
reduces space usage, but also reduces computation time. Our experimental
results showed that our algorithm is competitive with the state-of-the-art
tools under not so repetitive scenarios, and that it greatly reduces the
computational requirements when the input becomes more repetitive, standing
out as the most efficient tool to date (and to our knowledge) in this context.

## References

  * [1] Markus J. Bauer, Anthony J. Cox, and Giovanna Rosone.  Lightweight algorithms for constructing and inverting the BWT of string collections.  Theoretical Computer Science, 483:134–148, 2013.
  * [2] Hans-J Boehm, Russ Atkinson, and Michael Plass.  Ropes: An alternative to strings.  Software: Practice and Experience, 25(12):1315–1330, 1995.
  * [3] Paola Bonizzoni, Gianluca Della Vedova, Yuri Pirola, Marco Previtali, and Raffaella Rizzi.  Computing the multi–string BWT and LCP array in external memory.  Theoretical Computer Science, 862:42–58, 2021.
  * [4] Christina Boucher, Davide Cenzato, Zsuzsanna Lipták, Massimiliano Rossi, and Marinella Sciortino.  Computing the original eBWT faster, simpler, and with less memory.  In Proc. 28th International Symposium on String Processing and Information Retrieval (SPIRE), pages 129–142, 2021.
  * [5] Christina Boucher, Travis Gagie, Alan Kuhnle, Ben Langmead, Giovanni Manzini, and Taher Mun.  Prefix–free parsing for building big BWTs.  Algorithms for Molecular Biology, 14, 2019.  Article 13.
  * [6] Michael Burrows and David Wheeler.  A block sorting lossless data compression algorithm.  Technical Report 124, Digital Equipment Corporation, 1994.
  * [7] Diego Díaz-Domínguez and Gonzalo Navarro.  A grammar compressor for collections of reads with applications to the construction of the BWT.  In Proc. 31st Data Compression Conference (DCC), pages 83–92, 2021\.
  * [8] Lavinia Egidi, Felipe A. Louza, Giovanni Manzini, and Guilherme P. Telles.  External memory BWT and LCP computation for sequence collections with applications.  Algorithms for Molecular Biology, 14:6, 2019.
  * [9] Paolo Ferragina and Giovanni Manzini.  Opportunistic data structures with applications.  In Proc. 41st Annual Symposium on Foundations of Computer Science (FOCS), pages 390–398, 2000.
  * [10] Travis Gagie, Giovanni Manzini, and Jouni Sirén.  Wheeler graphs: A framework for BWT-based data structures.  Theoretical Computer Science, 698:67–78, 2017.
  * [11] Simon Gog, Timo Beller, Alistair Moffat, and Matthias Petri.  From Theory to Practice: Plug and Play with Succinct Data Structures.  In Proc. 13th International Symposium on Experimental Algorithms (SEA), pages 326–337, 2014.
  * [12] Richard Karp and Michael Rabin.  Efficient randomized pattern–matching algorithms.  IBM Journal of Research and Development, 31(2):249–260, 1987.
  * [13] Dominik Kempa.  Optimal construction of compressed indexes for highly repetitive texts.  In Proc. 30th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1344–1357, 2019.
  * [14] Dominik Kempa and Tomasz Kociumaka.  String Synchronizing Sets: Sublinear–Time BWT Construction and Optimal LCE Data Structure.  In Proc. 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC), page 756–767, 2019.
  * [15] Dominik Kempa and Tomasz Kociumaka.  Resolution of the Burrows–Wheeler transform conjecture.  In Proc. 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 1002–1013. IEEE, 2020.
  * [16] John C. Kieffer and En Hui Yang.  Grammar–based codes: a new class of universal lossless source codes.  IEEE Transactions on Information Theory, 46(3):737–754, 2000.
  * [17] Pang Ko and Srinivas Aluru.  Space efficient linear time construction of suffix arrays.  Journal of Discrete Algorithms, 3(2-4):143–156, 2005.
  * [18] Ben Langmead, Cole Trapnell, Mihai Pop, and Steven Salzberg.  Ultrafast and memory–efficient alignment of short DNA sequences to the human genome.  Genome Biology, 10(3), 2009.  Article R25.
  * [19] Jesper Larsson and Alistair Moffat.  Off–line dictionary–based compression.  Proceedings of the IEEE, 88(11):1722–1732, 2000.
  * [20] Heng Li.  Fast construction of fm-index for long sequence reads.  Bioinformatics, 30(22):3274–3275, 2014.
  * [21] Heng Li and Richard Durbin.  Fast and accurate long–read alignment with Burrows–Wheeler transform.  Bioinformatics, 26(5):589–595, 2010.
  * [22] Felipe A. Louza, Simon Gog, and Guilherme P. Telles.  Inducing enhanced suffix arrays for string collections.  Theoretical Computer Science, 678(1):22–39, 2017.
  * [23] Felipe A. Louza, Guilherme P. Telles, Simon Gog, Nicola Prezza, and Giovanna Rosone.  gsufsort: constructing suffix arrays, LCP arrays and BWTs for string collections.  Algorithms for Molecular Biology, 15, 2020.  Article 18.
  * [24] Robert Love.  Linux kernel development.  Pearson Education, 2010.
  * [25] Veli Mäkinen, Djamal Belazzougui, Fabio Cunial, and Alexandru Tomescu.  Genome-Scale Algorithm Design.  Camb. U. Press, 2015.
  * [26] Udi Manber and Gene Myers.  Suffix arrays: a new method for on–line string searches.  SIAM Journal on Computing, 22(5):935–948, 1993.
  * [27] Gonzalo Navarro.  Indexing Highly Repetitive String Collections, Part I: Repetitiveness Measures.  ACM Computing Surveys, 54(2), 2021.  Article 29.
  * [28] Gonzalo Navarro and Veli Mäkinen.  Compressed full-text indexes.  ACM Computing Surveys, 39:article 2, 2007.
  * [29] Ge Nong.  Practical linear-time O(1)-workspace suffix sorting for constant alphabets.  ACM Transactions on Information Systems, 31(3):1–15, 2013.
  * [30] Ge Nong, Sen Zhang, and Wai Hong Chan.  Linear suffix array construction by almost pure induced-sorting.  In Proc. 19th Data Compression Conference (DCC), pages 193–202, 2009.
  * [31] Daniel Saad Nogueira Nunes, Felipe A. Louza, Simon Gog, Mauricio Ayala-Rincón, and Gonzalo Navarro.  A grammar compression algorithm based on induced suffix sorting.  In Proc. 28th Data Compression Conference (DCC), pages 42–51, 2018\.
  * [32] Enno Ohlebusch.  Bioinformatics Algorithms: Sequence Analysis, Genome Rearrangements, and Phylogenetic Reconstruction.  Oldenbusch Verlag, 2013.
  * [33] Daisuke Okanohara and Kunihiko Sadakane.  Practical entropy-compressed rank/select dictionary.  In Proc. 9th Workshop on Algorithm Engineering and Experiments (ALENEX), pages 60–70. SIAM, 2007.
  * [34] Daisuke Okanohara and Kunihiko Sadakane.  A linear–time Burrows–Wheeler transform using induced sorting.  In Proc. 16th International Symposium on String Processing and Information Retrieval (SPIRE), pages 90–101, 2009.
  * [35] Zachary D. Stephens, Skylar Y. Lee, Faraz Faghri, Roy H. Campbell, Chengxiang Zhai, Miles J. Efron, Ravishankar Iyer, Michael C. Schatz, Saurabh Sinha, and Gene E. Robinson.  Big data: astronomical or genomical?  PLoS Biology, 13(7):e1002195, 2015.