analytical_paths <- list(
token_normalization = c("none", "lemmatization", "stemming"),
stopword_removal = c(TRUE, FALSE),
trimming = c(TRUE, FALSE),
alternative_model = c(TRUE, FALSE),
k_setting = c(1, 2, 3), #K original, alt1, alt2
iteration_setting = c(1, 2, 3) #iter original, alt1, alt2
)
settings_df <- expand.grid(analytical_paths, stringsAsFactors = FALSE)
nrow(settings_df)
#> [1] 2163 Multiverse analysis for topic modeling: A walkthrough of the procedure
Note
This is not a general tutorial on multiverse analysis. There are several guides on this already available. 1
This guide serves as an explanation to the logic of applying multiverse analysis for topic modeling used in the current study. We did not use any special multiverse analysis library to implement our analysis. Instead, we developed it almost from scratch in R, except using the functional programming library purrr and the hashing function rlang::hash(). This guide on one hand functions as the “annotated source code” to explain our code. On another hand, it helps researchers to understand multiverse analysis from first principles.
In this guide, we will use Tvinnereim and Fløttum (2015) as the example. Please note that unlike in the study, we do not consider “stochasticity” as an analytical path. If you would like to do that, please rerun the analysis. 2
3.1 Generate settings for all universes
Suppose all reasonable analytical paths are available (and have been preregistered), one can generate the Cartesian product of all analytical paths using expand.grid().
The output of expand.grid() is data frame. For easier manipulation, the data frame is converted to a list of lists called settings using purrr::transpose().
settings_df |> head()
#> token_normalization stopword_removal trimming alternative_model k_setting
#> 1 none TRUE TRUE TRUE 1
#> 2 lemmatization TRUE TRUE TRUE 1
#> 3 stemming TRUE TRUE TRUE 1
#> 4 none FALSE TRUE TRUE 1
#> 5 lemmatization FALSE TRUE TRUE 1
#> 6 stemming FALSE TRUE TRUE 1
#> iteration_setting
#> 1 1
#> 2 1
#> 3 1
#> 4 1
#> 5 1
#> 6 1settings <- settings_df |> purrr::transpose() settings is a list containing each combination of settings. settings and each singular item from settings can be subset in the following two ways. The subsetting of an item of settings is important because the first three slots from a set of settings are solely for text preprocessing.
# first two items
settings[1:2]
#> [[1]]
#> [[1]]$token_normalization
#> [1] "none"
#>
#> [[1]]$stopword_removal
#> [1] TRUE
#>
#> [[1]]$trimming
#> [1] TRUE
#>
#> [[1]]$alternative_model
#> [1] TRUE
#>
#> [[1]]$k_setting
#> [1] 1
#>
#> [[1]]$iteration_setting
#> [1] 1
#>
#>
#> [[2]]
#> [[2]]$token_normalization
#> [1] "lemmatization"
#>
#> [[2]]$stopword_removal
#> [1] TRUE
#>
#> [[2]]$trimming
#> [1] TRUE
#>
#> [[2]]$alternative_model
#> [1] TRUE
#>
#> [[2]]$k_setting
#> [1] 1
#>
#> [[2]]$iteration_setting
#> [1] 1# first item, first three slots
settings[[1]][1:3]
#> $token_normalization
#> [1] "none"
#>
#> $stopword_removal
#> [1] TRUE
#>
#> $trimming
#> [1] TRUE3.2 Generate all preprocessed document-term matrices
There are actually 12 ways to process the documents.
settings |> purrr::map(\(x) x[1:3]) |>
unique() |>
length()
#> [1] 12All these 12 document-term matrics will be generated.
preprocessing_settings <- settings |>
purrr::map(\(x) x[1:3]) |>
unique()3.2.1 Reading data
library(here)
library(readr)
library(quanteda)
library(udpipe)
library(dplyr)## the data is corrupted due to file encoding issues
input <- iconv(
readLines(here("rawdata", "ncp-stm-data.csv")),
from = "latin1",
to = "UTF-8"
) |>
paste(collapse = "\n") |>
readr::read_delim(delim = ";", show_col_types = FALSE) |>
na.omit()
original_corpus <- corpus(input, text_field = "openanswer")
original_corpus
#> Corpus consisting of 2,115 documents and 7 docvars.
#> text1 :
#> "Å være miljø vennlig "
#>
#> text2 :
#> "å så?har altid vært klimaendring "
#>
#> text3 :
#> "Å stange hodet mot veggen "
#>
#> text4 :
#> "å nei ikke igjen - gjør noe med det - NÅ"
#>
#> text5 :
#> "Ønsker ikke at klima skal endres "
#>
#> text6 :
#> "økte utfordringer med infrastrukturen i forhold til ekstremv..."
#>
#> [ reached max_ndoc ... 2,109 more documents ]3.2.2 Lemmatization
Lemmatization is done with udpipe.
norwegian_model <- udpipe_load_model(
file = here::here("rawdata/norwegian-bokmaal-ud-2.1-20180111.udpipe")
)
parsed_content <- udpipe_annotate(norwegian_model, original_corpus)
parsed_content_df <- as.data.frame(parsed_content)
parsed_content_df |>
group_by(doc_id) |>
summarize(content = paste(lemma, collapse = " ")) |>
mutate(doc_id = stringr::str_extract(doc_id, "[0-9]+")) |>
mutate(doc_id = as.numeric(doc_id)) |>
arrange(doc_id) -> lemma_df
lemma_corpus <- corpus(lemma_df$content)
docvars(lemma_corpus) <- docvars(original_corpus)
docnames(lemma_corpus) <- docnames(original_corpus)3.2.3 Token normalization
To realize all tokenization scenarios, three different ways of tokenization of the corpus were conducted: original (no lemmatization or stemming), lemmatization, and stemming. Please note that some code (e.g., the so-called “pretokenization”) was from the original code by Tvinnereim and Fløttum (2015), which converts some Norwegian terms and typos (e.g., globaloppvarming to global oppvarming, tempratur to temperatur).
3.2.3.1 Original (no lemmatization or stemming)
original_toks <- tokens(
original_corpus,
remove_punct = TRUE,
remove_numbers = TRUE,
include_docvars = TRUE
)3.2.3.2 Lemmatization
## the lemmatizer left $ before puntuation. Need to remove it explicitly
lemma_toks <- tokens(
lemma_corpus,
remove_punct = TRUE,
remove_numbers = TRUE,
include_docvars = TRUE
) |>
tokens_remove("$", valuetype = "fixed")3.2.3.3 Stemming
## Keep the original so-called "pre stemming" (very error prone, but we respect the original authors)
do_pre_stemming <- function(oa) {
for (i in 1:length(oa)) {
oa[[i]] <- gsub("frem", "fram", oa[[i]])
oa[[i]] <- gsub("forurensing", "forurens", oa[[i]])
oa[[i]] <- gsub("forurensning", "forurens", oa[[i]])
oa[[i]] <- gsub("smelting", "smelt", oa[[i]])
oa[[i]] <- gsub("altid", "alltid", oa[[i]])
oa[[i]] <- gsub("arktisk", "arktis", oa[[i]])
oa[[i]] <- gsub("bekymring", "bekymr", oa[[i]])
oa[[i]] <- gsub("bekymringsfullt", "bekymr", oa[[i]])
oa[[i]] <- gsub("betydning", "bety", oa[[i]])
oa[[i]] <- gsub("betyr", "bety", oa[[i]])
oa[[i]] <- gsub("død", "dø", oa[[i]])
oa[[i]] <- gsub("dør", "dø", oa[[i]])
# oa[[i]] <- gsub("endring", "endr", oa[[i]])
oa[[i]] <- gsub("enkelt", "enkel", oa[[i]])
oa[[i]] <- gsub("ekstremt", "ekstrem", oa[[i]])
oa[[i]] <- gsub("extrem", "ekstrem", oa[[i]])
oa[[i]] <- gsub("fleir", "fler", oa[[i]])
oa[[i]] <- gsub("flomm", "flom", oa[[i]])
oa[[i]] <- gsub("flomr", "flom", oa[[i]])
oa[[i]] <- gsub("forandring", "forandr", oa[[i]])
oa[[i]] <- gsub("fosilt", "fossil", oa[[i]])
oa[[i]] <- gsub("fossilt", "fossil", oa[[i]])
oa[[i]] <- gsub("fremtid", "framtid", oa[[i]])
oa[[i]] <- gsub("globaloppvarming", "global oppvarming", oa[[i]])
oa[[i]] <- gsub("godt", "god", oa[[i]])
oa[[i]] <- gsub("høgar", "høy", oa[[i]])
oa[[i]] <- gsub("høyer", "høy", oa[[i]])
oa[[i]] <- gsub("høyt", "høy", oa[[i]])
oa[[i]] <- gsub("konsekvens", "konsekv", oa[[i]])
oa[[i]] <- gsub("langt", "lang", oa[[i]])
oa[[i]] <- gsub("laver", "lav", oa[[i]])
oa[[i]] <- gsub("lavt", "lav", oa[[i]])
oa[[i]] <- gsub("meir", "mer", oa[[i]])
# oa[[i]] <- gsub("menneskeskapt", "menneskeskap", oa[[i]])
# oa[[i]] <- gsub("menneske", "mennesk", oa[[i]])
oa[[i]] <- gsub("overdrevent", "overdriv", oa[[i]])
oa[[i]] <- gsub("overdrev", "overdriv", oa[[i]])
oa[[i]] <- gsub("oson", "ozon", oa[[i]])
oa[[i]] <- gsub("ozonlag", "ozon", oa[[i]])
oa[[i]] <- gsub("politikern", "politiker", oa[[i]])
oa[[i]] <- gsub("reell", "reel", oa[[i]])
oa[[i]] <- gsub("reelt", "reel", oa[[i]])
oa[[i]] <- gsub("somr", "sommer", oa[[i]]) # exception: shorter to longer
oa[[i]] <- gsub("teknologisk", "teknologi", oa[[i]])
oa[[i]] <- gsub("temperaturendr", "temperaturforandr", oa[[i]])
oa[[i]] <- gsub("temperaturøkning", "temperaturstigning", oa[[i]])
oa[[i]] <- gsub("tempratur", "temperatur", oa[[i]]) # mis-spelling
oa[[i]] <- gsub("usikker", "usikk", oa[[i]])
oa[[i]] <- gsub("ustabilt", "ustabil", oa[[i]])
oa[[i]] <- gsub("utrydning", "utrydd", oa[[i]])
oa[[i]] <- gsub("utslepp", "utslipp", oa[[i]])
oa[[i]] <- gsub("uver", "uvær", oa[[i]])
oa[[i]] <- gsub("varmar", "varm", oa[[i]])
oa[[i]] <- gsub("varmerevåter", "varm våt", oa[[i]])
oa[[i]] <- gsub("varmer", "varm", oa[[i]])
oa[[i]] <- gsub("varmt", "varm", oa[[i]])
oa[[i]] <- gsub("vatn", "vann", oa[[i]])
oa[[i]] <- gsub("viktiger", "vikt", oa[[i]])
oa[[i]] <- gsub("viktigst", "vikt", oa[[i]])
oa[[i]] <- gsub("vinter", "vint", oa[[i]])
oa[[i]] <- gsub("vintr", "vint", oa[[i]])
oa[[i]] <- gsub("ødelagt", "ødel", oa[[i]])
oa[[i]] <- gsub("ødelegg", "ødel", oa[[i]])
oa[[i]] <- gsub("økend", "øke", oa[[i]])
oa[[i]] <- gsub("øker", "øke", oa[[i]])
oa[[i]] <- gsub("øket", "øke", oa[[i]])
oa[[i]] <- gsub("økning", "øke", oa[[i]])
oa[[i]] <- gsub("økt", "øke", oa[[i]])
oa[[i]] <- gsub("", "", oa[[i]])
oa[[i]] <- gsub("", "", oa[[i]])
}
return(oa)
}
prestemmed_toks <- do_pre_stemming(input$openanswer) |>
tokens(remove_punct = TRUE, remove_numbers = TRUE)
docvars(prestemmed_toks) <- docvars(original_toks)
docnames(prestemmed_toks) <- docnames(original_toks)3.2.3.4 Join the three together
current_tokens_list <- list()
current_tokens_list[["normal"]] <- original_toks
current_tokens_list[["lemmatized"]] <- lemma_toks
current_tokens_list[["prestemmed"]] <- prestemmed_toks3.2.4 Generate document-term matrices
process_tokens <- function(setting, current_tokens_list) {
if (setting$token_normalization == "lemmatization") {
current_tokens <- current_tokens_list[["lemmatized"]]
} else if (setting$token_normalization == "none") {
current_tokens <- current_tokens_list[["normal"]]
} else {
current_tokens <- current_tokens_list[["prestemmed"]]
}
if (setting$stopword_removal) {
current_tokens <- current_tokens |>
tokens_remove(
stopwords("norwegian"),
case_insensitive = TRUE,
padding = FALSE,
)
}
if (setting$token_normalization == "stemming") {
current_tokens <- tokens_wordstem(
current_tokens,
language = "norwegian"
)
}
current_dfm <- dfm(current_tokens)
if (setting$trimming) {
current_dfm <- dfm_trim(
current_dfm,
min_docfreq = 6,
docfreq_type = "count"
)
}
return(current_dfm)
}
dfm_list <- purrr::map(preprocessing_settings,
process_tokens,
current_tokens_list = current_tokens_list)In order to speed up the look up, each slot in dfm_list is named with the hash of the original preprocessing_setting.
names(dfm_list) <- purrr::map_chr(preprocessing_settings,
rlang::hash)So that one can quickly access the corresponding document-term matrix from a setting, e.g.,
dfm_list[[rlang::hash(settings[[100]][1:3])]]
#> Document-feature matrix of: 2,115 documents, 498 features (98.70% sparse) and 7 docvars.
#> features
#> docs å være miljø så har vært klimaendring mot ikke igjen
#> text1 1 1 1 0 0 0 0 0 0 0
#> text2 1 0 0 1 1 1 1 0 0 0
#> text3 1 0 0 0 0 0 0 1 0 0
#> text4 1 0 0 0 0 0 0 0 1 1
#> text5 0 0 0 0 0 0 0 0 1 0
#> text6 0 0 0 0 0 0 0 0 0 0
#> [ reached max_ndoc ... 2,109 more documents, reached max_nfeat ... 488 more features ]3.3 Train topic models
library(keyATM)
library(quanteda)
library(stm)The function train_model is built that accepts one setting and dfm_list and do the topic modeling for each setting in settings.
#' This function unifies the theta so that the output theta always
#' has the same nrow as current_dfm
#' The raison d'être is that stm discards empty rows silently
#' But doesn't retain the rownames
unify_theta <- function(mod, trimmed_dfm, current_dfm) {
theta <- mod$theta
rownames(theta) <- docnames(trimmed_dfm)
excluded_docs <- setdiff(docnames(current_dfm), docnames(trimmed_dfm))
k <- ncol(theta)
fake_theta <- matrix(
rep(1 / k, length(excluded_docs) * k),
nrow = length(excluded_docs),
ncol = k
)
rownames(fake_theta) <- excluded_docs
final_theta <- rbind(theta, fake_theta)
final_theta <- final_theta[
match(rownames(current_dfm), rownames(final_theta)),
]
return(final_theta)
}
train_model <- function(setting, dfm_list) {
current_dfm <- dfm_list[[rlang::hash(setting[1:3])]]
# iteration settings for original algorithm (STM) and
# alternative (Weighted LDA)
original_iter <- c(100, round(100 * 0.8), round(100 * 1.2))
alternative_iter <- c(2000, round(2000 * 0.8), round(2000 * 1.2))
## translate `setting` to actual topic model parameters `current`
k <- c(4 , 3, 5)[setting$k_setting]
if (!setting$alternative) {
iter <- original_iter[setting$iteration_setting]
} else {
iter <- alternative_iter[setting$iteration_setting]
}
random_seed <- sample(-65535:65535, 1)
output <- list()
output$random_seed <- random_seed
output$setting <- setting
set.seed(random_seed)
## STM drops rows silently, we do it here explicitly
rowsum_priv <- apply(current_dfm, 1, sum)
trimmed_dfm <- current_dfm[rowsum_priv != 0, ]
if (!setting$alternative_model) {
output$mod <- stm(
trimmed_dfm,
K = k,
init.type = "LDA", # See note above
max.em.its = iter,
prevalence = ~ concern + humanmade + efficacy + edu3 + gender + age,
verbose = FALSE,
data = trimmed_dfm@docvars
)
} else {
keyATM_docs <- keyATM_read(texts = trimmed_dfm)
output$mod <- weightedLDA(
docs = keyATM_docs,
number_of_topics = k,
model = "covariates",
model_settings = list(
covariates_data = trimmed_dfm@docvars,
covariates_formula = ~ concern +
humanmade +
efficacy +
edu3 +
gender +
age
),
options = list(
iterations = iter,
verbose = FALSE
)
)
}
output$theta <- unify_theta(output$mod, trimmed_dfm, current_dfm)
output$docvars <- trimmed_dfm@docvars
return(output)
}Then, the topic modeling for all settings is executed. This step roughly takes six hours.
models <- purrr::map(settings, train_model, dfm_list = dfm_list, .progress = TRUE)
names(models) <- purrr::map_chr(settings,
rlang::hash)3.4 Anchor
This setting represents the original setting used in Tvinnereim and Fløttum (2015). It is named anchor_setting.
anchor_setting <- list(token_normalization = "stemming",
stopword_removal = TRUE,
trimming = TRUE,
alternative_model = FALSE,
k_setting = 1,
iteration_setting = 1)This anchor setting is used to extract the model trained with this setting.
anchor_mod <- models[[rlang::hash(anchor_setting)]]And then the anchor model is used to perform exactly the same analysis that Tvinnereim and Fløttum (2015) did, namely, to find out the most distinguish topic cluster that is associated with age. And then extract the vector of topic proportion as anchor_theta (\(G_{:,A}\) in the article).
est <- stm::estimateEffect(
~ concern + edu3 + gender + age,
stmobj = anchor_mod$mod,
metadata = anchor_mod$docvars
)
max_topic_index <- which.max(
purrr::map_dbl(est$parameters, \(x) {
mean(purrr::map_dbl(x, \(y) y$est["age"]))
})
)
anchor_theta <- anchor_mod$theta[, max_topic_index]Suppose there is another \(\theta\) from another topic model.
another_mod <- models[[1]]
theta <- another_mod$thetatheta (\(\theta\)) is a matrix with \(k\) columns (total number of topics), that has topic proportions of all topics from a topic model solution resulting from a different set of settings. The idea of anchoring is to find out the column vector in \(\theta\) that has the maximum correlation with anchor_theta (See equation 2 in the article), which most likely represents the original author’s solution. So that one can use that column vector to represent the setting used for training this model.
cor_coefs <- vapply(
seq_len(ncol(theta)),
FUN = function(x) {
cor(anchor_theta, theta[, x], method = "spearman")
},
FUN.VALUE = numeric(1)
)
cor_coefs
#> [1] 0.70292271 0.05794823 -0.56911497 0.24747101Therefore, the column vector that can be used to represent this model is: theta[,which.max(cor_coefs)].
3.5 Generate estimates
With all the topic models trained, one can then generate estimates of the estimand. As stated in the paper, the estimand of interest is “the modeled difference in \(\theta_{:,k}\) across age groups.” To make it easier to do this for all trained topic models, the analyses are wrapped into functions.
#' return the column index in theta, which the column vector
#' has the highest spearman's correlation with anchor_theta
#' by default (return_rho = FALSE), it returns the anchor index
#' otherwise, it returns the maximun rho
find_anchor <- function(anchor_theta, theta, return_rho = FALSE) {
stopifnot(length(anchor_theta) == nrow(theta))
cor_coefs <- vapply(
seq_len(ncol(theta)),
FUN = function(x) {
cor(anchor_theta, theta[, x], method = "spearman")
},
FUN.VALUE = numeric(1)
)
if (return_rho) {
return(max(cor_coefs))
}
return(which.max(cor_coefs))
}
.get_keyatm_strata_topic_func <- function(current_mod) {
keyATM::by_strata_DocTopic(
current_mod$mod,
by_var = "age",
by_values = c(3, 4),
labels = c(3, 4)
)
}
.get_stm_estimate_func <- function(current_mod) {
stm::estimateEffect(
~ concern + humanmade + efficacy + edu3 + gender + age,
stmobj = current_mod$mod,
metadata = current_mod$docvars
) |>
plot(
covariate = "age",
model = current_mod$mod,
method = "difference",
cov.value1 = 4,
cov.value2 = 3,
omit.plot = TRUE
)
}
get_effect_size_mod <- function(
current_mod,
anchor_theta,
args,
.get_keyatm_strata_topic_func,
.get_stm_estimate_func
) {
setting <- current_mod$setting
set.seed(current_mod$random_seed)
k <- ncol(current_mod$theta)
if (setting$alternative_model) {
strata_topic <- .get_keyatm_strata_topic_func(current_mod)
theta1 <- strata_topic$theta[[1]]
theta2 <- strata_topic$theta[[2]]
theta_diff <- theta2[, seq_len(k)] - theta1[, seq_len(k)]
theta_diff_quantile <- apply(theta_diff, 2, quantile, c(0.025, 0.975))
theta_diff_mean <- apply(theta_diff, 2, mean)
anchor_index <- find_anchor(anchor_theta, current_mod$theta)
max_rho <- find_anchor(
anchor_theta,
current_mod$theta,
return_rho = TRUE
)
output <- data.frame(
Estimate = theta_diff_mean[anchor_index],
Q2.5 = theta_diff_quantile[1, anchor_index],
Q97.5 = theta_diff_quantile[2, anchor_index],
max_rho = max_rho
)
rownames(output) <- NULL
} else {
res <- .get_stm_estimate_func(current_mod)
anchor_index <- find_anchor(anchor_theta, current_mod$theta)
max_rho <- find_anchor(
anchor_theta,
current_mod$theta,
return_rho = TRUE
)
output <- data.frame(
Estimate = as.vector(res$means)[anchor_index],
Q2.5 = res$cis[[anchor_index]][1],
Q97.5 = res$cis[[anchor_index]][2],
max_rho = max_rho
)
colnames(output) <- c("Estimate", "Q2.5", "Q97.5", "max_rho")
rownames(output) <- NULL
}
output <- round(output, 6)
estimate <- cbind(as.data.frame(setting), output)
theta <- current_mod$theta[, anchor_index]
return(list(estimate = estimate, theta = theta))
}And then execute it via purrr::map.
estimates <- purrr::map(
models,
get_effect_size_mod,
anchor_theta = anchor_theta,
args = args,
.get_keyatm_strata_topic_func = .get_keyatm_strata_topic_func,
.get_stm_estimate_func = .get_stm_estimate_func)The output is joined as one data frame.
estimates_df <- estimates |>
purrr::map("estimate") |>
purrr::list_rbind()
head(estimates_df)
#> token_normalization stopword_removal trimming alternative_model k_setting
#> 1 none TRUE TRUE TRUE 1
#> 2 lemmatization TRUE TRUE TRUE 1
#> 3 stemming TRUE TRUE TRUE 1
#> 4 none FALSE TRUE TRUE 1
#> 5 lemmatization FALSE TRUE TRUE 1
#> 6 stemming FALSE TRUE TRUE 1
#> iteration_setting Estimate Q2.5 Q97.5 max_rho
#> 1 1 0.048773 0.028413 0.069289 0.702923
#> 2 1 0.055998 0.018252 0.082921 0.729928
#> 3 1 0.096883 0.074124 0.123890 0.656112
#> 4 1 -0.038074 -0.075956 0.002195 0.625493
#> 5 1 0.065408 0.034572 0.091677 0.684588
#> 6 1 0.077994 0.058560 0.104631 0.6766253.6 Visualize
For visualization, a customized visualization function is used. See lib/plot_spec_curve.R, which is based on the visualization function in the R package specr by Masur and Scharkow (2020)).
source(here::here("lib/plot_spec_curve.R"))tmmv.plot_spec_curve(estimates_df)
3.7 References
Masur, Philipp K., and Michael Scharkow. 2020. Specr: Conducting and Visualizing Specification Curve Analyses (Version 1.0.0). https://CRAN.R-project.org/package=specr.
Pipal, Christian, Hyunjin Song, and Hajo G. Boomgaarden. 2022. “If You Have Choices, Why Not Choose (and Share) All of Them? A Multiverse Approach to Understanding News Engagement on Social Media.” Digital Journalism, March, 1–21. https://doi.org/10.1080/21670811.2022.2036623.
Tvinnereim, Endre, and Kjersti Fløttum. 2015. “Explaining Topic Prevalence in Answers to Open-Ended Survey Questions about Climate Change.” Nature Climate Change 5 (8): 744–47. https://doi.org/10.1038/nclimate2663.
For example: the tutorial by Pipal et al. (2022) : https://cpipal.github.io/multiverse-tutorial/multiverse_tutorial.html↩︎
In comparison to the original code, we also remove many I/O operations, which are for optimization.↩︎