Quantitative predictions of neoadjuvant chemotherapy effects in breast cancer

by individual patient data assimililation

P. Castorina (a,b,∗), D.Carco’(a), C.Colarossi(a), M.Mare(a,e), L.Memeo(a)

,

M.Pace(b,c,d), I.Puliafito(a), D.Giuffrida(a)

a Dipartimento di Oncologia Sperimentale,Istituto Oncologico del Mediterraneo
, 95029 Viagrande, Italy

b INFN, Sezione di Catania, I-95123 Catania, Italy.

c Scuola di Specializzazione in Fisica Medica, Universita’ di Catania, Italy

d Centro Siciliano Fisica Nucleare e Struttura della Materia, Catania, Italy

e Dipartimento di Scienze Biomediche,odontoiatriche

e delle immagini morfologiche e funzionali, Universita’ di Messina,Italy

(*) corresponding author<EMAIL_ADDRESS>

Abstract

Neoadjuvant chemotherapy has been used for breast cancer aiming at downgrading
before surgery. In this article we propose a new quantitative analysis of the
effects of the neoadjuvant therapy to obtain numerical, personalized,
predictions on the shrinkage of the tumor size after the drug doses, by data
assimilation of the individual patient. The algorithm has been validated by a
sample of 37 patients with histological diagnosis of locally advanced primary
breast carcinoma. The biopsy specimen, the initial tumor size and its
reduction after each treatment were known for all patients. We find that: a)
the measure of tumor size at the diagnosis and after the first dose permits to
predict the size reduction for the follow up; b) the results are in agreement
with our data sample, within 10-20 $\%$, for about 90$\%$ of the patients. The
quantitative indications suggest the best time for surgery. The analysis is
patient oriented, weakly model dependent and can be applied to other cancer
phenotypes.

##  I Introduction

In the era of personalized oncology, mathematical models are a useful tool for
a better understanding of the clinical effects of therapy.

The expected individual response to tumor therapies is generally based on a
set of indices, defined with large quantitative variability. For example, for
neoadjuvant chemotherapy for locally advanced breast cancer, aiming at
downgrading before surgery, one usually considers the subtypes classification
according to the expression of hormone receptors, estrogen (ER) and
progesteron(PR), of the human epidermal growth factor receptor 2 (HER2), the
proliferation index ki67, the inizial tumor size and cellularity. Indeed,
these clinical informations may have a prognosis value similar to that of
multigene prognostic score cuzik .

The tumor progression during neoadjuvant chemotherapy claude ; fukuda ,
described by previous ( and others) predictive factors, gives direct
informations on the response to the therapy. By those analyses one gets semi-
quantitative results following the standard classification: tumor size (median
and range), T stage, Node stage, JACC stage, Lymphovascular invasion and other
parameters.

A complementary strategy could be obtained by more quantitative informations ,
based on numerical approaches which, by single patient data assimilation,
enhance the level of reliability of forecasts on the individual response.

Here we discuss an algorithm which, starting from the measure of the tumor
size ( radius ) at the diagnosis and after the first dose, is able to predict,
essentially without free parameters, the shrinkage of the tumor in the
sequence of treatments. The proposed method is , by itself, patient oriented
since the first size reduction and the initial cellularity take into account
the specific initial condition.

The numerical predictions agree , within $10-20\%$, for more than $90\%$ of
the observed data of our sample of 37 patients.

Figure 1: Predictions to data ratio on the tumor radius after the second dose.
The green (red) line indicates the $1\pm 0.1$ ($1\pm 0.2$) interval.

##  II Mathematical formulation of the diagnostic algorithm

The breast tumor growth is described by the Gompertz law gompertz ; norton1 ,
solution of the differential equation

$({1}/{N})dN/dt=-kln({N}/{N_{\infty}}),$ (1)

where $N$ is the cell number at time $t$, $k$ is a constant and $N_{\infty}$
is the maximum number of cells ($N_{\infty}=3.1*10^{12}$ , according to ref.
norton1 ).

The modification of the specific growth rate due to chemotherapy, during the
time interval of a single treatment, is obtained by introducing a function
$c(t)$ in the previous equation norton2 ; norton3 ; norton4 ; castorina , i.e.

$({1}/{N})dN/dt=-kln({N}/{N_{\infty}})-c(t),$ (2)

where $c(t)$ has a negligible value after the interval, $\tau$, between two
timeline doses ($\tau=3$ weeks). In other terms, chemoterapy effects start,
periodically, at the beginning of each drug dose but almost completely decline
after $\tau=3$ weeks and, therefore, the function $c(t)$ has a discontinuity
on the days of treatment. By solving the previous eq.(2) ( see appendix A) for
homogeneous, spherical symmetric configurations, the size reduction after $n$
doses is given by

$ln({R_{n+1}}/{R_{0}})=ln({R_{1}}/{R_{0}})\Sigma_{m=0}^{n}exp{(-mk\tau)},$ (3)

where $R_{n+1}$ is the tumor radius after $n+1$ doses and, for each patient,
the constant $k$ is determined by the initial cellularity (the second term in
the growth law in eq.(1) is the fraction of duplicating tumor cells).

In the final result, (eq.3), the function $c(t)$ does not explicitely appear:
its contribution is hidden in the (measured) size after the first dose
$R_{1}(\tau)$. In this respect, the approach is independent on the model
describing the chemotherapy effects.

##  III Validation, Results and Discussion

###  III.1 Validation: Patients and Therapy

Patients This is a retrospective single centre study. Thirty-seven women, aged
36-78 years, with histologically proven operable breast cancer were evaluated.
All tumours were tested for estrogen receptor (ER), progesterone receptor
(PgR), HER 2 and ki 67 . Thirty-six patients showed positivity for ER (range
$2-90\%$) and PgR (range $2-90\%$), HER2 $3+$ was present in $5/37$ patients.
Ki 67 was variable from $5$ to $30\%$. Median diameter of tumour , defined by
imaging, was $43,5$ mm (range $21-72$ mm) .Four patients had clinical
positivity for axillary nodes. Pregnant women were excluded. ECOG-PS of all
patients was 0 or 1. All patients had adequate haematological, renal and
haepatic function. All patients had a normal left ventricular ejection
fraction (LVEF > $50\%$). Treatment Neoadjuvant chemotherapy corresponds to
the use of a systemic treatment applied before locoregional treatment (surgery
and /or radiotherapy) in order to obtain a more frequent conservating surgery,
downgrading the tumour size. Major drugs used for breast cancer patients
included anthracyclines and taxanes iv1  . Patients evaluated in our study
received a median of five cycles (range 4-6) of every -3-week(q3w) ET
(epirubicin 80 mg /m2 , paclitaxel 175 mg/m2 ) iv2 ; iv3 . Seven patients
having HER2 $3+$ received integrate treatment with trastuzumab 6 mg/kg ( 8
mg/kg as loading dose). At the first follow up , after one chemotherapy
administration , all patients had a tumour diameter reduction variable from
$10$ to $70\%$. At the second follow up, after second chemotherapy
administration , all patients showed a further diameter reduction included
between $10$ and $30\%$. At the third follow up, $14/37$ patients continued to
respond to treatment while the others showed a stable disease . At the fourth
follow up , only one patient showed a futher tumour diameter reduction , the
others continued to have a stabilization of disease and this was persisting at
the remaining follow up iv4 .

###  III.2 Results and Discussion

The estimate of the tumor shrinkage in the dose sequence follows immediately
from eq.(3) and from the determination of $R_{1}(\tau)$. In Fig.1 and Fig.2
the numerical results are compared with data for the second and the third
treatment for all patients. The radius measurements had a $2-3\%$ statistical
error and the error propagation has been taken into account.

For the second dose, the ratio between predictions and data is within the
prudential interval $1\pm 0.1$ ($1$ indicates a perfect agreement) for 31
patients of the entire sample ($84\%$) and the agreement is within $1\pm 0.2$
(see Fig.1) for $34/37$ ( $92\%$).

For the third dose, in $29/37$ and $32/37$ cases the agreement is within the
fiducial $1\pm 0.1$, $1\pm 0.2$ intervals (see Fig.2) respectively.

The results in Fig.2 have been obtained by assuming an almost constant tumor
size for stable disease. On the other hand, one can ask if the diagnostic
algorithm can give quantitative indications on the stability of the disease,
i.e. if there is only a small reduction of the size after the treatment.

If a further reduction less than $10\%$ defines the stable disease condition,
by applying the proposed algorithm, one gets that after the third treatment,
$17/34$ patients continued to respond to the therapy (i.e. the tumor size
decreases more than $10\%$) to be compared with the clinical result of our
sample, $14/37$ cases ( see Fig.3).

In Fig.4 analogous results are reported for the fourth dose, giving $4/37$
patients still responding to the treatment ($1/37$ is the clinical result).

Figure 2: Predictions to data ratio on the tumor radius after the third dose.
The green (red) line indicates the $1\pm 0.1$ ($1\pm 0.2$) interval.

A homogeneous and spherical tumor is assumed in the previous sections, however
these constraints can be easily removed and the diagnostic algorithm can be
improved in many directions, taking into account, for example, different
geometrical tumor shapes or inhomogeneities.

The proposed approach can be applied in the clinical practice as follows.
Initially the tumor size, $R_{0}$ and the cellularity are evaluated by usual
methods and, after the first dose, one measures the radius ( $R_{1}$). By
those input data, one then estimates the shrinkage of the tumor size due to
the following treatments, i.e the values $R_{2},R_{3},..$, according to
eq.(3). If the tumor shrinkage observed after the second dose ( and before the
third one) turns out to be in agreement with or larger than the estimated
results , $R_{2}$, then the predictions should be considered clinically
reliable. In particular, a forecast of stable disease ( defined by a reduction
of the tumor radius less than $10\%$) suggests to stop the neo-adjuvant
therapy and to proceed with surgery: there is a clear signal that other drug
doses are not effective to further reduce the size.

A computational code can be easily implemented.

##  IV Conclusions

The good agreement between predictions and data for the treatment sequence
suggests that the proposed method is a reliable starting point for a more
quantitave description of neo-adjuvant chemotherapy effects and for an optimal
management of patients, permitting to avoid unnecessary treatments and
reducing economic costs. It should be further clarified that it has to be
considered as a complementary tool to the standard chemotherapy response
evaluation criteria in solid tumors iv1 ; viale ; colleoni ; CPSEG ; recist1 ;
recist2 ; choi  and that it does not give any information on the overall
survival probability, but a quantitative evaluation of the tumor size
depletion.

Figure 3: Predictions of stabe disease after the third dose. The red lines
indicate the $1\pm 0.1$ interval. Figure 4: Predictions of stabe disease after
the fourth dose. The red lines indicate the $1\pm 0.1$ interval.

Acknowledgments

The authors thank the Fondazione IOM. The work is partially supported by the
project "DiOncoGen Diagnostica Innovativa", Codice CUP: G89J18000700007,
Azione 1.1.5 del PO FESR SICILIA 2014/2020

## References

  * (1) J.Cuzik et al., Prognostic value of a combined ER,PgR,Ki67,HER2 immunohistochemical (IHC4) score and comparison with GHI recurrence score-results from TransATAC, Cancer Res 2009,69:503s
  * (2) Abigail S. Claude et al.,Predictors of tumor progression during neoadjuvant chemotherapy in breast cancer, J. Clin. Oncol. 28:1821-1828 (2010)
  * (3) I.Fukuda et al., Pattern of tumor shrinkage during neoadjuvant chemotherapy is associated with prognosis in low grade luminal early breast cancer, Radiology 2018,286:49-57
  * (4) Pathak et al. Neoadjuvant chemotherapy regimens in treatment of breast cancer: a systematic review and network meta-analysis protocolSystematic Reviews (2018) 7:8
  * (5) von Minckwitz G. Docetaxel, Anthracycline combinations for breast cancer treatment. Expert Opin Pharmacother. 2007;8(4):485–95.
  * (6) Zaheed M, Wilcken N, Willson ML, O Connell DL, Goodwin A., Sequencing of anthracyclines and taxanes in neoadjuvant and adjuvant therapy for early breast cancer. Cochrane Database Syst Rev. 2019 Feb 18;2:CD012873
  * (7) Caitlin Murphy et al. Tailored NEOadjuvant epirubicin, cyclophosphamide and Nanoparticle Albumin-Bound paclitaxel for breast cancer: The phase II NEONAB trial—Clinical outcomes and molecular determinants of response. PLOS ONE | https://doi.org/10.1371/journal.pone.0210891 February 14, 2019
  * (8) B. Gompertz , On the nature of the function expressive of t he law of human mortality and a new mode of determining life contingencies”, Phil. Trans. R. Soc. 115,513 (1825)
  * (9) L. Norton ,A. Gompertzian model of human breast cancer gro wth” Cancer Res 48,7067-7071 (1988)
  * (10) Norton L, Simon R., Tumor size, sensitivity to therapy and the design of treatment protocols. Cancer Treat Rep 1976;61:1307–1
  * (11) Norton L., Conceptual and practical implications of breast cancer geometry: toward a more effective, less toxic therapy. Oncologist 2005;10:370–8
  * (12) Norton L, Simon R. ,The Norton-Simon hypothe- sis: designing more effective and less toxic chemo- therapeutic regimens. Nat Clin Pract Oncol 2006;3: 406–7
  * (13) P. Castorina et al., Tumor Growth Instability and Its Implications for Chemotherapy, Cancer Res 2009;69:8507-8515.
  * (14) Viale G, Regan MM, Maiorano E, et al. Prognostic and predictive value of centrally reviewed expression of estrogen and progesterone receptors in a randomized trial comparing letrozole and tamoxifen adjuvant ther- apy for postmenopausal early breas cancer: BIG 1-98. J Clin Oncol 2007;25:3846–52
  * (15) Colleoni M. et al. Expression of ER, PgR, HER1, HER2, and response: a study of preop- erative chemotherapy. Ann Oncol 2008;19:465–72
  * (16) Mittendorf E.A. et al. , The Neo-Bioscore Update for Staging Breast Cancer Treated With Neoadjuvant Chemotherapy-Incorporation of Prognostic Biologic Factors Into Staging After Treatment. JAMA Oncology, 2(7), 929. doi:10.1001/jamaoncol.2015.6478
  * (17) Therasse P., Arbuck S.G., Eisenhauer E.A., et al. New guidelines to evaluate the response to treatment in solid tumors. European Orga- nization for Research and Treatment of Cancer, National Cancer Institute of the United States, National Cancer Institute of Canada. J Natl Cancer Inst. 2000;92:205-216.
  * (18) Eisenhauer E.A., Therasse P., Bogaerts J., et al. New response evalua- tion criteria in solid tumors: revised RECIST guideline (version 1.1). Eur J Cancer. 2009;45:228-247.
  * (19) Choi H., Charnsangavej C., Faria S.C., at al. Correlation of computed tomography and positron emission tomography in patients with metastatic gastrointestinal stromal tumor treated at a single institu- tion with imatinib mesylate: proposal of new computed tomography response criteria. J Clin Oncol. 2007;25:1753-1759.

##  V Appendix A

The Gompertz law is solution of the equation gompertz ,

$\frac{1}{N}\frac{dN}{dt}=-kln(\frac{N}{N_{\infty}}),$ (4)

which describes the macroscopic growth of a cancer cell population, without
the effect of chemotherapy norton1 . $N(t)$ is the cell number at time t, $k$
is a constant and $N_{\infty}$ is the maximum number of cells ( carrying
capacity).

The drug treatment modifies the previous equation by inroducing the function
$c(t)$ which decribes the chemotherapy effects during a single treatment, i.e.
it is different from zero in the interval, $\tau$, between two timeline doses
($\tau=3$ weeks). The general solution of the equation turns out to be

$y(t)=y(0)e^{-k(t-t_{0})}-e^{-k(t-t_{0})}\int_{t_{0}}^{t}c(t^{{}^{\prime}})e^{+k(t^{\prime}-t_{0})}dt^{{}^{\prime}}$
(5)

where $y(t)=ln(N(t)/N_{\infty})$, and $y(0)$ is the value at the initial time
of the treatment $t_{0}$. By defining the initial time $t_{0}=0$ and

$I(\tau)=e^{-k\tau}\int_{0}^{\tau}c(t^{{}^{\prime}})e^{+kt^{\prime}}$ (6)

the result of the first dose is given by

$y_{1}(\tau)=y(0)e^{-k\tau}-I(\tau)$ (7)

Since the function $c(t)$ is different from zero only for time interval $\tau$
after the beginning of any single dose, the initial condition for the second
treatment is $y_{1}(\tau)$ and at the end ($t=2\tau$) one gets

$y_{2}(2\tau)=y_{1}(\tau)e^{-k\tau}-I(\tau)$ (8)

By using the previous eq.(5) one obtains

$y_{2}(\tau)-y(0)=[y_{1}(\tau)-y(0)][e^{-k\tau}+1]$ (9)

and, by iteration, at the end of $n+1$ treatments, one immediately gets

$y_{n+1}[(n+1)\tau]-y(0)=[y_{1}(\tau)-y(0)]\Sigma_{m=0}^{n}e^{-mk\tau}$ (10)

The previous equation written in term of tumor volume, for uniform density, is

$ln(\frac{V_{n+1}}{V_{0}})=ln(\frac{V_{1}}{V_{0}})\Sigma_{m=0}^{n}e^{-mk\tau}$
(11)

where $V_{n+1}$ is the tumor volume after $n+1$ doses. If one considers
spherical symmetry, the previous equation can be written as

$ln(\frac{R_{n+1}}{R_{0}})=ln(\frac{R_{1}}{R_{0}})\Sigma_{m=0}^{n}e^{-mk\tau}$
(12)

where $R_{n+1}$ is the tumor radius after $n+1$ doses.


# Efficient implementation of sets and multisets in R using hash tables

Giacomo Ceoldo, Ernst C. Wit

###### Abstract

The package hset for the R language contains an implementation of a S4 class
for sets and multisets of numbers. The implementation, based on the hash table
data structure from the package hash (Brown, 2019), allows for quick
operations when the set is a dynamic object. An important example is when a
set or a multiset is part of the state of a Markov chain in which in each
iteration various elements are moved in and out of the set.

##  1 Introduction

Sets are the most basic and fundamental containers of objects in mathematics.
According to set theory (almost) all objects in mathematics are, or can be
described as, sets. Some objects have additional mathematical and
computational structure, such as multisets, lists, vectors, stacks, etc. Sets
and multisets have been less developed in programming languages than other
objects, such as vectors and lists. The reason is that the latter ones are
used very frequently in algorithms, so almost all programming languages have a
built-in implementation of them. These implementations do not reflect the
mathematical derivation of these objects from set theory, as their structure
allows the use of much more computationally efficient implementations and
algorithms. Sets and multisets as basic containers are nevertheless very
important, especially for discrete probabilistic and statistical models. The
aim of the paper is to provide an efficient implementation for algorithms that
use sets and multisets as containers. Our implementation is based on the hash
package (Brown, 2019). It is efficient because the hash table data structure
allows search, insertion and deletion of one element in the table in constant
time.

The mathematical definition of sets and multisets is given in Section 2, where
our implementation and its semantic is also discussed. In Section 3, relations
and operations between sets and multisets are defined mathematically, and
their implementation is described. The performance of our implementation is
discussed in Section 4. In Section 5 an application of our implementation of
sets and multisets as states of a Markov chain is provided.

##  2 Sets and multisets

Sets, multisets and some other constructions/containers derived from them, are
introduced mathematically in Section 2.1. The emphasis is on how containers
differ on how much structure is “imposed” on them. In sets elements are either
in or out of them, in multisets it is also relevant how frequently an element
is contained, for sequences the order of the elements is also important. The R
package hset is introduced in Section 2.2, where it is also discussed which
objects can be included as elements, and how they are stored in a data
structure based on the package hash. The semantic of hset objects is described
in Section 2.3, with some other functions that are used to control these
objects, that depend on the chosen semantic.

###  2.1 Mathematical definition

_Sets_ are defined as collections of objects called _elements_ , or _members_.
_Set theory_ can be used as foundation of mathematics. The existence of the
_empty set_ $\emptyset=\\{\\}$ is postulated, $\emptyset$ is the only element
such that

$\\{a,\emptyset\\}=\\{a\\},\ \forall a,$ (1)

where equality between sets will be formally defined in Section 3.1. All
elements of a set are considered to be sets themselves, and (almost) all
objects in mathematics are constructed as sets.

An important example is the set $\mathbb{N}$ of _natural numbers_ that can be
defined recursively as

$0=\emptyset,\ 1=\\{0\\}=\\{\emptyset\\},\
2=\\{0,1\\}=\\{\emptyset,\\{\emptyset\\}\\},\ ...$ (2)

Another important example is the _Cartesian product_ of two sets, that is the
set of ordered pairs with first element from $X_{1}$ and second element from
$X_{2}$, defined as

$X_{1}\times X_{2}=\\{X_{1},\\{X_{2}\\}\\},$ (3)

which contains elements $(a_{1},a_{2})$, for $\emptyset\neq a_{1}\in X_{1}$
and $\emptyset\neq a_{2}\in X_{2}$. Note that $\emptyset\times
X=X\times\emptyset=\\{\emptyset\\}$ for all $X$, and that
$(a_{1},a_{2})\neq(a_{2},a_{1})$ for $a_{1}\neq a_{2}$ (the Cartesian product
is not _commutative_).

Order and multiplicity of the elements of a set is not defined, that is

$\\{a,b\\}=\\{b,a\\}=\\{a,a,b\\},\ \forall a\neq b.$ (4)

_Multisets_ are defined as collections of elements with _multiplicities_ ,
that are non-negative numbers. We will always assume that the multiplicities
of the elements are finite. The order of the elements is not defined, but
multisets of elements with different multiplicities are different:

$\\{a[m],b[n]\\}\neq\\{a[n],b[m]\\},\ \forall n\neq m,\ \forall\emptyset\neq
a\neq b\neq\emptyset$ (5)

If the multiplicity of an element is $0$, then the element is not contained in
the multiset:

$\\{a[m],b[0]\\}=\\{a[m]\\},\ \forall a,b.$ (6)

A set can be considered equivalent to a multiset of the same elements, all
with multiplicity $1$:

$X=\\{a,b\\}\cong\\{a[1],b[1]\\}=Y(X),\ \forall a,b,$ (7)

The injective function $X\mapsto Y(X)$ converts a set to the “equivalent”
multiset, but when applied to a multiset, it is the identity function:
$Y=Y(Y)$. The surjective function $Y\mapsto X(Y)$ maps the multiset $Y$ into
its _support_ , that is the set $X(Y)$ of its elements, if $X$ is a set
$X=X(X)$. The size of a multiset (or of a set), is the number of its elements,
that is

$\mathrm{size}(Y)=|Y|=|X(Y)|\in\mathbb{N}.$ (8)

The cardinality of a multiset $Y=\\{a_{i}[m_{i}]\\}_{i}$ is the sum of its
multiplicities, that is

$\mathrm{card}(Y)=||Y||={\textstyle\sum}_{i}m_{i}\in[0,\infty),$ (9)

while for a set $X$, $\mathrm{size}(X)=\mathrm{card}(X)$.

Other constructions that will be used later on are based on the Cartesian
product, that is used to define _powers_ of the set $X$ as

$X^{k}=X\times X\times...\times X,$ (10)

for $k\in\mathbb{N}$, where $X$ is repeated $k$ times in the right hand side,
$X^{0}=\emptyset$ and $X^{1}=X$. Powers of $X$ are used themselves to define
finite dimensional _sequences_ (or _strings_) of elements in $X$ as

$A=(a_{1},...,a_{l})\in X^{*}=\biguplus_{k\geq 0}X^{k},$ (11)

so a sequence of length $l$ with elements in $X$ is an element of $X^{l}$, and
the set of sequences $X^{*}$, that is the _disjoint union_ of all powers,
contains all sequences of all finite possible lengths, including the _empty
sequence_ $()$. Sequences have even more structure than multisets, as they are
ordered, meaning that $(a_{1},a_{2},...)\neq(a_{2},a_{1},...)$ for $a_{1}\neq
a_{2}$.

Finite dimensional sequences are introduced because they are one of the most
used objects in programming, so most languages have efficient built-in
implementations of sequences with finite length. However these implementations
are often inefficient when a sequence is modified locally by operations of
inclusion or removal of elements from the list.

In R finite dimensional sequences are implemented as objects of type vector,
with sub-types atomic and list (Wickham, 2019, Chapter 3). In the next
sections, our R implementation of sets and multisets is discussed in detail.

###  2.2 Computational implementation

In mathematics, elements of sets are usually considered sets themselves. Then
a formal definition of elements is redundant, once sets are defined. On the
other hand, when sets are viewed as computational objects, a definition of
elements is required because the elements are objects with, possibly various,
data types stored in memory.

Sets and multisets are implemented in the _R package_ hset, as objects of _S4
class_ "hset". These objects are containers of elements, that are either
numbers, or sets of numbers. A more formal definition of objects that are
valid elements will soon be given. An object of class "hset" contains two
_slots_. The main one, called @htable, is an hash table from the package hash
Brown (2019). The second slot, called @info, is of class "environment", which
contains a Boolean value that distinguish sets from multisets. The reason why
the second slot is a “trivial” (with one object) environment, rather than an
object of class "logical", is because environments and (logical) vectors have
a different semantic, then it would be difficult to reason about “sameness” of
objects (sets or multisets). The constructor for objects of class "hset" will
be described at the end of this section, the semantic of our implementation
will be discussed in detail in the next section.

The set S of possible sets that can be stored is recursively defined as

$\begin{split}&X=\\{a_{1},a_{2},...\\}\in\texttt{S},\\\
&a_{i}\in\texttt{S}\uplus\texttt{N},\end{split}$ (12)

so the element $a_{i}$ can be either a set, or a value in N that is the set of
numeric vectors of length $1$, without the values Inf, -Inf, NaN, NA,
NA_integer_ and NA_real_, that are excluded. The symbol $\uplus$ in
$\texttt{S}\uplus\texttt{N}$ is used to denote the disjoint union between the
sets S and N. The inclusion relation $\in$ between an element, and a set or
multiset, will be formally defined in Section 3.1.

The set M of possible multisets that can be stored is recursively defined as

$\begin{split}&Y=\\{a_{1}[m_{1}],a_{2}[m_{2}],...\\}\in\texttt{M},\\\
&a_{i}\in\texttt{S}\uplus\texttt{N},\\\
&m_{i}\in\texttt{N}_{\texttt{+}},\end{split}$ (13)

where S and N are defined as above, and
$\texttt{N}_{\texttt{+}}\subset\texttt{N}$ is the subset of values of N that
are strictly positive. Note that in our current implementation multisets can
not be elements of a set or a multiset, and that recursion occurs in the
definition of M because it occurs in the definition of S. The numeric datatype
includes integer and double as subtypes, so the elements can be also of these
two types. Vectors of type numeric of length 0, together with the NULL object,
are considered equivalent to the empty set, so they are not included in an
object of class "hset" (even though formally an empty set is included in every
set). Instead, numeric values of length at least $2$ and list values of every
length are converted to elements of S, that is to sets, before being included
as elements.

The package sets (Meyer and Hornik, 2009) uses a different approach, the sets
in this library, that can be of classes set (sets), gset (_generalized_ sets),
cset (_customizable_ sets), can contain elements of every type. Two elements
can be considered the same only if they have the same class, so for example
the sets $\\{\texttt{2},\texttt{2L}\\}$ and $\\{\texttt{2}\\}$, are not equal
(the former has two elements), as 2 and 2L have class "numeric" and "integer"
respectively. As a result, although our implementation is more limited because
we only take into account sets and multisets of numbers, it is still somewhat
closer to the mathematical definition in which 2L and 2 represent the same
number.

The _hash table_ implemented in hash is a data structure that contains _key-
value pairs_. The _keys_ are different objects of type character, they are
unique labels of pointers to the _values_ , that can be objects of every type
except NULL. The advantage of using an hash table to implement sets and
multisets is that various operations that are: adding and removing a key-value
pair from the table, checking to see if a key is present in the table, and
returning the value associated with a key; only require on average a constant
number of elementary operations.

There is an injection $k:\texttt{S}\uplus\texttt{N}\to\texttt{C}$, where C is
the set of character vectors of length 1, that is used as set of keys to label
uniquely each possible element of a set or a multiset. For example, the
element

$a_{i}=\\{\texttt{-1},\texttt{1},\texttt{1L},\texttt{\\{\\}},\texttt{2},\texttt{11},\texttt{\\{2,\\{3\\}\\}}\\},$
(14)

is mapped to

$k_{i}=k(a_{i})=\texttt{"\\{-1,\\{\\},\\{2,\\{3\\}\\},1,11,2\\}"},$ (15)

where the ”sub-elements” 1 and 1L are both mapped to "1", and the components
of the character $k_{i}$ are in _lexicographic order_ , which has the
mathematical properties of a _total order_. Note that the order is important
to guarantee that $k$ is an injection. For sets, all values of the hash table
are equal to the empty character "", that is $v(k_{i})=\texttt{""}$ for all
$i$, whereas for multisets $v(k_{i})\in\texttt{N}_{\texttt{+}}$.

Sets and multisets are created with the _constructor_ hset with three
arguments that are members, multiplicities and generalized. The first two
arguments are NULL by default (empty set), the last one is FALSE by default.
If the second input is not NULL, generalized is set to TRUE if it was not the
case. Then it is checked that members and multiplicities are of the correct
type, that they are coherent with themselves, and then they are included in
the hash table. The function is.hset, with input x, returns TRUE when x is of
class "hset", and FALSE otherwise. The function as.hset, with input x, return
x itself if it is of class "hset", otherwise it applies the constructor hset,
with members equal to x, multiplicities and generalized as default, so that
the function creates a set with elements taken from x.

Size and cardinality defined in equations (8) and (9) are returned by the
functions size.support and cardinality. A vector containing the labels of the
elements is returned by the function members, while the vector of
multiplicities of the elements are returned by the function multiplicities.
The only input of these four functions is an object of class "hset". The
components of the vectors obtained by the last two functions are coherent, so
that the $i$-th value of the vector is the multiplicity of the $i$-th element.
If the function multiplicities is used on a set, a vector with all values
equal to $1$ is returned.

###  2.3 Semantic of hset

In R, objects can be accessed and modified with _reference_ , or with _value
semantic_ (the two alternatives are described in Appendix A). In R objects of
most classes have value semantic, but environments and hash tables from the
package hash have reference semantics. An object of class "hset" contains a
hash object in the first slot, and an environment object in the second. When a
set is copied “directly”, both slots of the copied object, refer to the same
ones of the original object. The functions clone.of.hset and refer.to.hset are
used to copy an hset with value and reference semantic, respectively.

The function is.generalized, with logical returned value, is used to
distinguish sets and multisets. The map as.generalized transforms a set to a
multiset by converting the values of all elements of the hash table to 1L,
while as.not.generalized transforms a multisets to a set by converting all
values to "". These functions implement $X\mapsto Y(X)$ and $Y\mapsto X(Y)$
that were defined mathematically in Section 2.1. The sets are modified
locally, i.e., with reference semantic, so if a multiset is transformed to a
set the information about the multiplicities of the multiset that is passed as
input is lost. The functions clone.of.hset and refer.to.hset can also be used
to convert sets to multisets, and vice versa. They have as second argument the
(empty, or logical) value called generalized, that is NULL by default, but it
can be used to convert a multiset to a set, or a set to a multiset, when it is
copied. Applying refer.to.hset with second argument equal to TRUE
(respectively FALSE), is equivalent to applying the function as.generalized
(respectively as.not.generalized). Whereas the application of clone.of.hset
with second argument equal to TRUE or FALSE, creates a new hset with the same
support as the original one, but with the multiplicities that are converted to
1L or "" in the two cases, and the hset that is passed as input of
clone.of.hset is not modified.

In the next section relations and operations between sets and multisets will
be described mathematically and computationally. In this section we describe
how the chosen semantic can affect the computation of an operation. Relations
are encoded as functions with Boolean codomain, so if two components that can
be elements, sets or multisets, are in relation, the function returns TRUE,
otherwise FALSE. As sets are not modified when checking whether two components
are in relation, the semantic is irrelevant. Conversely for operations between
sets or multisets, even though the result of the operation is not changed by
the semantic used, one operand will be modified when reference semantic is
used, while a new hset with the result of the operation is created. The
operands do not change when value semantic is used.

All operations are computed with the function hset.operation.numeric if at
least one of the operands is a multiset, otherwise hset.operation.logical is
used. These two functions have the same signature in which the output is an
object called new.hset of class "hset". The first input hset1 is the first
operand, ... contains all other operands (for operations with multiple arity),
the arguments operation (function) and identity.is.universe (logical value)
completely specify how the operands are combined. The last input semantic,
that can be equal to "refer" (default) or "value", specifies the semantic. The
difference between the numeric and the logical operation is in how the
multiplicities are handled. In the latter case we define the multiplicities of
a set by the bijection that maps NULL to FALSE and "" to TRUE, where the
Boolean outcomes are used to evaluate the operation. The evaluation is stored
using the inverse function, as setting an element to NULL in an hash object is
equivalent to removing the key-value pair from the hash table, or to not doing
anything if the pair is not present. If the numeric operation is used, the
multiplicities that are used in the operation are obtained by the surjective
function of type
$\texttt{NULL}\uplus\texttt{""}\uplus\texttt{N}_{\texttt{+}}\to\texttt{N}_{\texttt{+=}}$,
s.t. $\texttt{NULL}\mapsto\texttt{0}$, $\texttt{""}\mapsto\texttt{1}$, and
$\texttt{m}\mapsto\texttt{m}$ for $\texttt{m}\in\texttt{N}_{\texttt{+}}$. The
non-negative values are then used as operands, and the result is stored in the
hash table with the bijection
$\texttt{N}_{\texttt{+=}}\to\texttt{NULL}\uplus\texttt{N}_{\texttt{+}}$, s.t.
$\texttt{0}\mapsto\texttt{NULL}$, and $\texttt{m}\mapsto\texttt{m}$ for all
$\texttt{m}\in\texttt{N}_{\texttt{+}}$. Note that the output of a numeric
operation is always a multiset, so "" is never stored.

The function create.new.hset is used in hset.operation.numeric and
hset.operation.logical, to create the object new.hset that will store the
result. When reference or value semantics are used, refer.to.hset or
clone.of.hset are used respectively, inside create.new.hset, with argument
hset1. Therefore, new.hset and hset1 will refer to the same object in memory
with reference semantic, so that when the result is computed, will be stored
both in new.hset, and in hset1. Whereas with value semantic, new.hset will be
a reference to an object in memory that is a clone of hset1, so the latter
will not change when the result is computed. The reference semantic is used by
default for its computational advantages. In particular when the identity
element of the operation is the empty set, that is when identity.is.universe
is FALSE, computing the result does not require a complete scan through the
elements of hset1.

The computational complexity of an operation between multisets $Y_{1}$,
$Y_{2}$, $...$, $Y_{a}$, where $\emptyset$ is the identity element of the
operation, is $O(|Y_{2}|+...+|Y_{a}|)$ with reference semantic, and
$O(|Y_{1}|+|Y_{2}|+...+|Y_{a}|)$ with value semantic. The advantage is
significant when $|Y_{1}|\gg\max_{j\neq 1}|Y_{j}|$. Note that difference of
$O(|Y_{1}|)$ operations between the two semantics, is due to the necessity of
copying the first operand. However, if the first operand is a set, while some
of the others are multisets, there is no difference between the two semantics,
because the result of the operation is a multiset, so $O(|Y_{1}|)$ operations
are required to convert the first operand to a multiset, even when reference
semantic is used.

##  3 Sets algebra

Relations and operations involving hset objects are described in Sections 3.1
and 3.2, respectively. Relations of different types are described by their
signature, definition, and implementation as functions with Boolean codomain.
Operations are also described in the same way, but the functions that compute
them return sets or multisets.

###  3.1 Relations

#### Inclusion of elements.

The _inclusion_ relation between an element and a set is defined
mathematically as

$\in:(\texttt{N}\uplus\texttt{S})\times\texttt{S},\ \ \ \ a\in X\ \iff\
X=\\{a,...\\},$ (16)

meaning that the element $a\in(\texttt{N}\uplus\texttt{S})$ and the set
$X\in\texttt{S}$ are related by $\in$ if and only if $a$ is an element of $X$.
Note that $\emptyset\in X$ for all sets $X$, and that the symbol $:$ is used
in the signature of the relation to avoid the notation
$\in\in(\texttt{N}\uplus\texttt{S})\times\texttt{S}$. This relation is
extended trivially to multisets as

$\in:(\texttt{N}\uplus\texttt{S})\times\texttt{M},\ \ \ \ a\in Y\ \iff\
Y=\\{a[n],...\\},\ n\geq 1,$ (17)

that is if the multiplicity of $a\in(\texttt{N}\uplus\texttt{S})$ in
$Y\in\texttt{M}$ is at least $1$.

The straightforward generalization for multisets are relations between an
element with a given multiplicity, and a multiset. Three relations of this
type are defined as

$\in_{\sim}:((\texttt{N}\uplus\texttt{S})\times\texttt{N}_{\texttt{+}})\times\texttt{M},\
\ \ \ a[m]\in_{\sim}Y\ \iff\ Y=\\{a[n],...\\},\ m-n\sim 0,$ (18)

where $\sim$ can be $\leq$, $<$ and $=$, for the relations of, inclusion
$\in=\in_{\leq}$, _strict inclusion_ $\in_{<}$ and _exact inclusion_
$\in_{=}$, respectively. Intuitively, in the three cases, $a[m]$ and $Y$ are
related if $a$ in $Y$ has a multiplicity greater or equal, greater, and equal
to $m$ respectively. Instead of defining relations between an element with a
given multiplicity and a multiset, we could have equivalently defined the
family of relations between elements and multisets parametrized by
$m\in\texttt{N}_{\texttt{+}}$:

$\in_{\sim}^{m}:(\texttt{N}\uplus\texttt{S})\times\texttt{M},\ \ \ \
a\in_{\sim}^{m}Y\ \iff\ Y=\\{a[n],...\\},\ m-n\sim 0,$ (19)

so that, for all $a\in\texttt{N}\uplus\texttt{S}$, $Y\in\texttt{M}$, and
$m\in\texttt{N}_{\texttt{+}}$,

$a[m]\in_{\sim}Y\ \iff\ a\in_{\sim}^{m}Y.$ (20)

All relations defined above can be encoded as one function with signature

$((\texttt{N}\uplus\texttt{S})\uplus((\texttt{N}\uplus\texttt{S})\times\texttt{N}_{\texttt{+}}))\times(\texttt{S}\uplus\texttt{M})\times\\{\leq,<,=\\}\to\\{\texttt{TRUE},\texttt{FALSE}\\},$
(21)

where the first argument is either an element, or a pair between an element
and a multiplicity, the second argument is either a set or a multiset, whereas
the last argument specifies the type of relation. The function returns TRUE if
the first two arguments are in relation, of the type specified by the third
argument. In the package, a similar function, called inclusion.member, has
signature

$\texttt{C}\times(\texttt{S}\uplus\texttt{M})\times\texttt{N}_{\texttt{+}}\times\\{\leq,<,=\\}\to\\{\texttt{TRUE},\texttt{FALSE}\\}.$
(22)

The last two arguments, called multiplicity and type.relation with default
values 1 and $\leq$ respectively, are ignored when the second argument is a
set. The first argument, called member must be a vector of length 1, that is
converted to a character C inside the function, if this value is not a valid
element, as defined above, the function returns FALSE for every possible
choice of the last two arguments. Then, inclusion.member returns TRUE if and
only if the first two arguments are in relation specified by the last two
arguments. The binary operator %in% uses inclusion.member where the last two
arguments are set as default, so that it evaluates the relations (16) or (17),
depending on whether the second argument is a set or a multiset respectively.
If the first argument, i.e., the left operand of %in%, is a vector of
characters, the operand returns a vector of Booleans with the result of the
evaluated relation for each character.

#### Subsets and equalities.

Now relations in which both objects are sets or multisets are considered. The
_subset relation_ between sets $X_{1}$ and $X_{2}$ is

$\subseteq:\texttt{S}\times\texttt{S},\ \ \ \ X_{1}\subseteq X_{2}\iff(a\in
X_{1}\implies a\in X_{2}),$ (23)

and the _strict subset relation_ is

$\begin{split}\subset:\texttt{S}\times\texttt{S},\ \ \ \ X_{1}\subset
X_{2}&\iff(X_{1}\subseteq X_{2}\ \mathrm{and}\ \exists b\in X_{2}\
\mathrm{s.t.}\ b\notin X_{1}).\end{split}$ (24)

The _equality relation_ between sets is defined as

$\begin{split}=:\texttt{S}\times\texttt{S},\ \ \ \
X_{1}=X_{2}&\iff(X_{1}\subseteq X_{2}\ \mathrm{and}\ X_{2}\subseteq
X_{1}).\end{split}$ (25)

For finite dimensional sets, the last two relations can also be written as

$X_{1}\approx X_{2}\iff X_{1}\subseteq X_{2}\ \mathrm{and}\
|X_{1}|\sim|X_{2}|,$ (26)

where for the strict inclusion, $\approx$ and $\sim$ are replaced by $\subset$
and $<$ respectively, while for the equality relation $\approx$ and $\sim$ are
both replaced by $=$.

The relations above will be generalized in the case in which at least one
component of the relation is a multiset. For
$Y=\\{a_{i}[m_{i}]\\}\in\texttt{M}$, the multiplicities are written as a
function $v_{Y}:X(Y)\to\texttt{N}_{\texttt{+}}$, such that
$v_{Y}(a_{i})=m_{i}$. The domain of this function, that is the support of $Y$,
is extended to all well defined elements, as

$v_{Y}:\texttt{S}\uplus\texttt{N}\to\texttt{N}_{\texttt{+}}\cup\\{0\\},\ \
\text{s.t.}\ \ v_{Y}(a_{i})=\begin{cases}m_{i}&\text{if}\ a_{i}\in X(Y)\\\
0&\text{otherwise}\end{cases}.$ (27)

Note that the extension of this function is coherent with equation (6). The
relations above are generalized as

$\begin{split}&Y_{1}\subseteq Y_{2}\iff v_{Y_{1}}(a)\leq v_{Y_{2}}(a)\ \forall
a\in\texttt{S}\uplus\texttt{N},\\\ &Y_{1}\subset Y_{2}\iff Y_{1}\subseteq
Y_{2}\ \mathrm{and}\ \exists a\in\texttt{S}\uplus\texttt{N}\ \mathrm{s.t.}\
v_{Y_{1}}(a)<v_{Y_{2}}(a),\\\ &Y_{1}\sqsubseteq Y_{2}\iff Y_{1}\subseteq
Y_{2}\ \mathrm{and}\ \nexists a\in\texttt{S}\uplus\texttt{N}\ \mathrm{s.t.}\
0\neq v_{Y_{1}}(a)<v_{Y_{2}}(a),\\\ &Y_{1}\sqsubset Y_{2}\iff Y_{1}\sqsubseteq
Y_{2}\ \mathrm{and}\ \exists a\in\texttt{S}\uplus\texttt{N}\ \mathrm{s.t.}\
0=v_{Y_{1}}(a)<v_{Y_{2}}(a),\\\ &Y_{1}=Y_{2}\iff Y_{1}\subseteq Y_{2}\
\mathrm{and}\ Y_{2}\subseteq Y_{1}\iff v_{Y_{1}}(a)=v_{Y_{2}}(a)\ \forall
a\in\texttt{S}\uplus\texttt{N},\end{split}$ (28)

The signature of these relation is
$\approx:\texttt{S}\uplus\texttt{M}\times\texttt{S}\uplus\texttt{M}$, where
$\approx$ is one of the five relations above. However, only the definition for
multisets is given in (28), but this is not a problem, as equation (7) implies
that if at least one argument, say the first one, is a set, the relation
$X_{1}\approx Y_{2}$ is equivalent to $Y(X_{1})\approx Y_{2}$.

For finite dimensional multisets, the relations can be written as

$\begin{split}&Y_{1}\subseteq Y_{2}\iff v_{Y_{1}}(a)\leq v_{Y_{2}}(a)\ \forall
a\in X(Y_{1}),\\\ &Y_{1}\subset Y_{2}\iff Y_{1}\subseteq Y_{2}\ \mathrm{and}\
(\exists a\in X(Y_{1})\ v_{Y_{1}}(a)<v_{Y_{2}}(a)\ \mathrm{or}\
|Y_{1}|<|Y_{2}|),\\\ &Y_{1}\sqsubseteq Y_{2}\iff v_{Y_{1}}(a)=v_{Y_{2}}(a)\
\forall a\in X(Y_{1})\ \mathrm{and}\ |Y_{1}|\leq|Y_{2}|,\\\ &Y_{1}\sqsubset
Y_{2}\iff v_{Y_{1}}(a)=v_{Y_{2}}(a)\ \forall a\in X(Y_{1})\ \mathrm{and}\
|Y_{1}|<|Y_{2}|,\\\ &Y_{1}=Y_{2}\iff v_{Y_{1}}(a)=v_{Y_{2}}(a)\ \forall a\in
X(Y_{1})\ \mathrm{and}\ |Y_{1}|=|Y_{2}|.\end{split}$ (29)

Other definitions are available, but these ones are the most efficient
computationally, because it is not necessary to evaluates multiplicities in
$Y_{2}$ for elements not in $X(Y_{1})$. The function hset1.included.to.hset2
with signature

$(\texttt{S}\uplus\texttt{M})\times(\texttt{S}\uplus\texttt{M})\times\\{\texttt{TRUE},\texttt{FALSE}\\}\times\\{\texttt{TRUE},\texttt{FALSE}\\}\to\\{\texttt{TRUE},\texttt{FALSE}\\}$
(30)

where the four arguments are hset1, hset2, strictly and exactly, returns TRUE
if the first two arguments are in one of the relations defined above, except
for the equality, that is implemented with another function. If the first two
arguments are both sets, the fourth argument is ignored because in
$\texttt{S}\times\texttt{S}$ the relations $\subseteq$ and $\sqsubseteq$ are
equivalent, and so are $\subset$ and $\sqsubset$. The function iterates
through the members of hset1, computes the difference of multiplicities, if
this difference is negative FALSE is returned immediately, otherwise the
difference is accumulated. For $\sqsubseteq$ and $\sqsubset$ the accumulated
difference must be zero at the end of the loop, and the two relations are
distinguished by the condition on the supports. For $\subset$, either the
accumulated difference is strictly positive, or the support of the second set
is larger. Whereas for $\subseteq$, no other conditions are required after the
end of the loop. The evaluation of the equality relations is implemented in
the function hsets.are.equal, in which only the functions size.support,
members and multiplicities defined at the end of Section 2.1 are used to
evaluate the relation.

Some generic operators that call the function hset1.included.to.hset2 with
different combinations of the last two inputs are defined. For $\subseteq$, <=
and >= are used, where the latter is for the inverse relation $\supseteq$,
obtained by reflecting the arguments. For $\subset$, the generic operators are
< and >, for $\sqsubseteq$ they are %=<=% and %=>=%, and for $\sqsubset$, %=<%
and %=>% are used. The equality operator == calls the function hsets.are.equal
and != returns its negation.

###  3.2 Operations

When discussing relations the semantic of the implementation was never
mentioned, as the sets that were possibly be part of some relations were never
modified, and returned by the functions used to check whether two objects are
in relation. An _operation_ is a ternary relation between sets or multisets,
that is written as $X_{1}\approx X_{2}=X$ for a given $\approx$, however the
interest here is computing $X$ from the _operands_ $X_{1}$ and $X_{2}$.
Moreover all operations will be defined for general arity, that is for more
than two operands. The _universe set_ , denoted by $U$ is a set such that
$X\subseteq U$, for all $X\in\texttt{S}$, that is the set containing all
possible elements: $a\in U$ for all $a\in\texttt{S}\uplus\texttt{N}$. The
_universe multiset_ $U$ is the multiset such that $Y\subseteq U$, for all
$Y\in\texttt{M}$. Then, $U$ contains all elements in
$\texttt{S}\uplus\texttt{N}$, each with infinite multiplicity. Note that the
universe set and multiset are defined by the same symbol, as it will be clear
from the context which “universe” is considered.

The signature of all operations defined below is

$\approx:(\texttt{S}\uplus\texttt{M})\times\biguplus_{k\geq
0}(\texttt{S}\uplus\texttt{M})^{k}\to(\texttt{S}\uplus\texttt{M}),$ (31)

so that the operation is defined if there is at least one operand, the second
argument is written as a disjoint union of all possible tuple of hsets, so
that for a given $k$, there are $k+1$ operands. However, only binary
operations are defined, as the extension to general arities is
straightforward. The result of the operation is of class S if and only if all
operands are of class S, otherwise the result is of class M. In the latter
case, if an operand $X_{i}$ is a set, it is replaced by $Y_{i}=Y(X_{i})$.

The _intersection_ , _union_ , _sum_ , _difference_ and _symmetric difference_
between $X_{1}$ and $X_{2}$ are

$\begin{split}&X=X_{1}\cap X_{2}\iff(a\in X\implies a\in X_{1}\ \text{and}\
a\in X_{2}),\\\ &X=X_{1}\cup X_{2}\iff(a\in X\implies a\in X_{1}\ \text{or}\
a\in X_{2}),\\\ &X=X_{1}+X_{2}\iff X=X_{1}\cup X_{2},\\\ &X=X_{1}\setminus
X_{2}\iff(a\in X\implies a\in X_{1}\ \text{and}\ a\notin X_{2}),\\\
&X=X_{1}\triangle X_{2}\iff X=(X_{1}\setminus X_{2})\cup(X_{2}\setminus
X_{1}),\end{split}$ (32)

respectively. The identity element for the intersection is $X_{2}=U$, while
for all other operations is $X_{2}=\emptyset$. All operations, except for the
difference, are commutative and associative, $X_{1}\setminus X_{2}\setminus
X_{3}$ is defined to be equal to $(X_{1}\setminus X_{2})\setminus X_{3}$. Note
that the sum has been defined to be the same as the union, but these
operations will be different for multisets. The multiset version of the
operations above is

$\begin{split}&Y=Y_{1}\cap Y_{2}\iff\forall
a,v_{Y}(a)=\min(v_{Y_{1}}(a),v_{Y_{2}}(a)),\\\ &Y=Y_{1}\cup Y_{2}\iff\forall
a,v_{Y}(a)=\max(v_{Y_{1}}(a),v_{Y_{2}}(a)),\\\ &Y=Y_{1}+Y_{2}\iff\forall
a,v_{Y}(a)=v_{Y_{1}}(a)+v_{Y_{2}}(a),\\\ &Y=Y_{1}\setminus Y_{2}\iff\forall
a,v_{Y}(a)=\max(v_{Y_{1}}(a)-v_{Y_{2}}(a),0),\\\ &Y=Y_{1}\triangle
Y_{2}\iff\forall a,v_{Y}(a)=|v_{Y_{1}}(a)-v_{Y_{2}}(a)|.\end{split}$ (33)

As for the set version, the identity elements are $Y_{2}=U$ and
$Y_{2}=\emptyset$ for the intersection, and for all other operations
respectively. When generalizing $\setminus$ and $\triangle$ to multisets, some
properties that hold for sets are violated. For example, in general
$(Y_{1}\setminus Y_{2})\cap Y_{2}\neq\emptyset$, and $Y_{1}\triangle
Y_{2}\triangle Y_{3}\neq Y_{2}\triangle Y_{1}\triangle Y_{3}$.

The functions hset.operation.numeric and hset.operation.logical have signature

$\begin{split}&(\texttt{S}\uplus\texttt{M})\times(\texttt{S}\uplus\texttt{M})^{*}\times((\texttt{N}_{\texttt{+}}\uplus\texttt{0})^{*}\to\texttt{N}_{\texttt{+}}\uplus\texttt{0})\times\texttt{B}\times\\{\texttt{"refer"},\texttt{"value"}\\}\to\texttt{M},\\\
&\texttt{S}\times\texttt{S}^{*}\times(\texttt{B}^{*}\to\texttt{B})\times\texttt{B}\times\\{\texttt{"refer"},\texttt{"value"}\\}\to\texttt{S},\end{split}$
(34)

respectively, where $\texttt{Z}^{*}$ is computed as in equation (11), and
$\texttt{B}=\\{\texttt{FALSE},\texttt{TRUE}\\}$. The first two arguments,
called hset1 and ..., contain the operands. The third argument, called
operation, is a function that computes the updated multiplicity, the choice of
this function defines which of the operations above (intersection, union, sum,
difference, symmetric difference) is computed. The fourth argument called
identity.is.universe, with default FALSE, specifies whether iterating through
all elements of the first operand is needed. The last argument called
semantic, with default "refer", specifies whether the first operand is
modified, or whether its clone is modified. For the intersection, union, sum,
difference, and symmetric difference, with logical multiplicities, the
functions that are used in the third input are all, any, any, nimp (for “not
implies”) and niff (for “not if and only if”) respectively. Whereas, with
numeric multiplicities, the functions are min, max, sum, pdif (for “positive
difference”) and sdif (for “symmetric difference”) respectively. The functions
nimp, niff, pdif and sdif have been implemented by us, while the others are
primitive functions in R. In the intersection the fourth argument is TRUE ($U$
is the identity element of the operation), while in all other operations, the
fourth argument is FALSE ($\emptyset$ is the identity element).

The functions intersection, union, setsum, difference, symmdiff with signature

$(\texttt{S}\uplus\texttt{M})\times(\texttt{S}\uplus\texttt{M})^{*}\times\\{\texttt{"refer"},\texttt{"value"}\\}\to(\texttt{S}\uplus\texttt{M}),$
(35)

call hset.operation.numeric if at least one operand is a multiset, otherwise
hset.operation.logical is used, with third and fourth arguments set by the
operation. The infix operands for computing the binary intersection are %&%,
%&&% and %and%, for the union, the operands are %|%, %||% and %or%, for the
sum are %+% and %sum%, for the difference %-% and %!implies%, and for the
symmetric difference %--% and %xor%. In all these operands reference semantic
is used, for a binary operation with value semantic, all operands can be used,
but a $\sim$ is added before the last %, e.g., for the union with value
semantic, %|$\sim$% can be used.

##  4 Performance

Here we assess the performance of our implementation, by comparing its two
semantics between themselves, and with the package sets (Meyer and Hornik,
2009).

###  4.1 Relations

Computing a relation results to a Boolean output indicating whether the two
components are related. Therefore the sets are not modified and the semantic
is irrelevant. The only comparison is then between our implementation based on
an hash table, denoted here by hsets, and the implementation of the library
sets.

The comparison for the inclusion relation $a\in X$ from equation (16) between
the element $a$ and the set $X$, is in Figure 1. The _x_ axis is $|X|$, which
is the size of the set $X$, on the _y_ axis the time of evaluating the
relation ai %in% X, where ai is a vector of elements, and X have classes hset
and sets, for the black and the red dots, respectively. The shape of the dots
denotes for how many elements, that is the size of ai, the relation is
evaluated. In our implementation, the complexity of the operation does not
depend on the size of the set $|X|$, but it depends linearly on the size of
ai, that is the number of elements that we want to check whether they are
contained in the set. On the other hand, the complexity of the implementation
of sets, depends linearly on $|X|$, but it seems not to depend on the size of
ai, implying that the inclusion relations onto a set are parallellized in
sets.

Figure 1: Time complexity for evaluating the inclusion relation.

In our implementation, the subset and the equality relations between sets,
written succinctly as $X_{1}\approx X_{2}$, with $\approx$ in
$\\{\subseteq,\subset,=\\}$, have the same complexity evaluated above, that is
constant with respect to $|X_{2}|$, and linear with respect to $|X_{1}|$,
where $X_{2}=X$ and $X_{1}=\\{a_{i}\\}_{i}$. The same complexity can be found
if $X_{1}$ and $X_{2}$ are replaced by multisets $Y_{1}$ and $Y_{2}$, for all
cases of $\approx$. The reason is that, when evaluated using the formulas in
(29), the complexity of the relation grows linearly with respect to $|Y_{1}|$
because of the component $\forall a\in X(Y_{1})$, and the complexity is
constant with respect to $|Y_{2}|$ because $|Y_{1}|\approx|Y_{2}|$ is
evaluated in constant time, for various $\approx$.

###  4.2 Operations

The codomain of an operation is either a set or a multiset. In Figure 2 the
complexity of the defined set operations is plotted in the left column. The
_x_ and _y_ axes distinguish the size of the operands and the time for
evaluating one operation in seconds. The shape of the points denotes the size
of the second operands, their colours distinguish the semantics of our
implementation and the operations computed with objects from the package sets.
In Section 2.3, it has been explained that the reference semantic is helpful
for some operations when the first operand is large. In particular, the time
complexity of the operation is constant with respect to the size of the first
operand, when $\emptyset$ is the identity element of the operation, that is
for the union, difference and symmetric difference. The complexity cost of
using the value semantic can be seen by the fact that the complexity grows
linearly with the size of both operands, as the (large) first one has to be
copied at the beginning of the operation. Whereas with the reference semantic
the complexity of these operation is linear only with respect to the size of
the second operand, while being constant with respect to the size of the first
one. In the library sets, the complexity of the operations grows linearly with
respect to the size of the first operand, but it seems not to depend on the
size of the second one. We think that it actually depends on the size of the
largest operands. For the interaction however our implementation is much more
inefficient that sets, regardless of the semantic, although all
implementations have a linear complexity with respect to the size of the first
operand. The sum of two sets is defined to be equivalent to the union, so the
complexity of the implementations is not computed.

In the right column of Figure 2 the same comparison has been done for
operations between multisets. As in the previous case, the intersection is
much more efficient in the sets package, whereas for all other operations, our
implementation with reference semantic does not depend on the size of the
first operand, while with value semantic it does, as the first operand is
copied.

Figure 2: Time complexity of set (left) and multiset (right) operations for
different implementations, and sizes of the operands.

##  5 MCMC with state space of undirected graphs

In this section we describe an application of the hset package to _Markov
processes_ with set- or multiset-valued (discrete) _state space_. The example
considered here is a stochastic process with set-valued state space, and
multiset-valued _sufficient statistic_ for the distribution of the process.
The stochastic process is defined over the $n(n-1)/2$ dimensional state space
of $n$ dimensional _undirected simple graphs_. Each graph is represented
uniquely by the _edge set_ , which contains all information (ties) about the
graph. The state of the process is augmented with the _degree distribution_ of
the graph, that is a multiset and contains all relevant information about the
distribution of the stochastic process.

The network process evolves by _tie flips_ in which (typically few) non-edges
become edges and vice versa. The tie flips are often _local_ , e.g., when only
one tie can be flipped at given time, in general at each time point only a
number of ties that is small, in comparison with the size of the network, can
be flipped. If few edges of the graph are flipped, the sufficient statistic
does not have to be recomputed from the edge set, as only the degrees of the
vertices adjacent to the flipped edges must be updated. Therefore the
reference semantic in our implementation that is derived from the hash table
data structure, has the computational advantage of being able to modify the
state locally (in memory), without having to copy the state when its size
changes.

We consider here the _Beta Model_ (Holland and Leinhardt, 1981; Blitzstein and
Diaconis, 2011; Rinaldo et al., 2013), where vertex $i$ has its own parameter
$\beta_{i}$, and the tie $ij$ between vertices $i$ and $j$ is in the graph
with probability $p_{ij}=e^{\beta_{i}+\beta_{j}}/(1+e^{\beta_{i}+\beta_{j}}).$
The model is equivalent to a sample of $n(n-1)/2$ independent Bernoulli random
variables $X_{ij}$ with probabilities $p_{ij}$. The distribution of the
network can be written as

$P(X=x|\beta)=\prod_{1\leq i<j\leq
n}p_{ij}=\exp\Big{(}\sum_{i=1}^{n}d_{i}(x)\beta_{i}-\psi(\beta)\Big{)},$ (36)

where $d_{i}(x)$ is the degree of vertex $i$, and
$\psi(\beta)=\sum_{i<j}\log(1+e^{\beta_{i}+\beta_{j}})$ is the log-partition
function. Thus the _degree vector_ $d=(d_{1},...,d_{n})$ is a sufficient
statistic for the model.

We consider a Markov chain in which at each time point some (typically few, in
comparison with $n$) tie flips are proposed, and accepted with probability

$Q(\tilde{x}|x,\beta)=\min\Big{(}\exp\Big{(}{\sum}_{i\in\mathcal{I}}(\tilde{d}_{i}-d_{i})\beta_{i}\Big{)},1\Big{)},$
(37)

where $\mathcal{I}\subseteq\\{1,...,n\\}$ and $i\in I$ if exist vertex $j$
such that the tie $ij$ is included in the proposed tie flips (thus
$i\in\mathcal{I}\iff j\in\mathcal{I}$), $d_{i}=d_{i}(x)$ and
$\tilde{d}_{i}=d_{i}(\tilde{x})$ are the degrees of vertex $i$ before and
after the tie flips, respectively. If the proposed tie flip contains only the
tie $ij$, then $\mathcal{I}=\\{i,j\\}$ and $|\mathcal{I}|=2$, in general more
than one tie can be flipped, but the algorithm is useful computationally when
$|\mathcal{I}|\ll n$. Note that this process is a Markov chain with
_Metropolis_ updates, as we use a symmetric proposal distribution for the tie
flips (uniform distribution over $n(n-1)/2$). Therefore, the stationary
distribution of the Markov chain with acceptance probability
$Q(\tilde{x}|x,\beta)$ in equation (37) is $P(X=x|\beta)$ in equation (36).

Three Markov chains $(X_{t})_{t}$, $({X_{t}}^{-})_{t}$ and $({X_{t}}^{+})_{t}$
are considered, with same transition probability parametrized by $\beta$, but
with different starting points. The chains $(X_{t})_{t}$, $({X_{t}}^{-})_{t}$
and $({X_{t}}^{+})_{t}$ have _stationary_ , _sparse_ and _dense starting_
point, respectively. Therefore $X_{0}$, ${X_{0}}^{-}$ and ${X_{0}}^{+}$ have
degrees similar, lower and higher, respectively, than the expected degrees
computed from the stationary distribution of the Markov chain. The real
parameter is generated as $\beta\sim\mathrm{Norm}(-1_{n},I_{n})$. At each
iteration, a single tie flip is proposed, sampling
$h\sim\mathrm{Unif}(1,...,n(n-1)/2)$, and the tie $ij$ is computed from $h$.

With our implementation, the state of the chain with stationary starting point
is $(X_{t},Z_{t})$, that is composed by two objects of class "hset", $X_{t}$
is the network encoded as set of edges, and $Z_{t}$ is the degree distribution
of the network encoded as a multiset. In each iteration the hset of flips $F$
is sampled, in our case $|F|=1$ as only one tie is sampled. Then the set
$\mathcal{I}_{F}$ of vertices that are part of at least one proposed tie flip
is computed, and so are the proposed degrees $\tilde{d}_{i}$ for
$i\in\mathcal{I}_{F}$. If the flip is accepted, the state is updated as

$\begin{split}&X_{t+1}\xleftarrow{}X_{t}\ \triangle\ F,\\\
&Z_{t+1}\xleftarrow{}(Z_{t}-\\{d_{i}[m_{i}]:i\in\mathcal{I}_{F}\\})+\\{\tilde{d}_{i}[\tilde{m}_{i}]:i\in\mathcal{I}_{F}\\},\end{split}$
(38)

where $\tilde{m}_{i}$ and $m_{i}$ are the multiplicities of proposed
$\tilde{d}_{i}$ and current $d_{i}$ degrees respectively, for
$i\in\mathcal{I}_{F}$. The update uses the operations
$\triangle:\texttt{S}\times\texttt{S}\to\texttt{S}$,
$+:\texttt{M}\times\texttt{M}\to\texttt{M}$ and
$-:\texttt{M}\times\texttt{M}\to\texttt{M}$, that have been defined in Section
3.2. In Section 4.2 it has been shown that for these three operations (that
all have the empty set / multiset as the identity second operand), the
complexity is not affected by the sizes of $X_{t}$ and $Z_{t}$, but it depends
linearly on the sizes of $F$ and $\mathcal{I}_{F}$. The same algorithm is used
with $({X_{t}}^{+},{Z_{t}}^{+})$ and $({X_{t}}^{-},{Z_{t}}^{-})$, in which the
starting point is out of equilibrium .

The updates in equation (38) are coded with reference semantic as

$\begin{split}&\texttt{state\$edge.set \%xor\% hset(flips\$id)}\\\
&\texttt{state\$degree.frequencies \%-\% hset(names(table.old.degrees),}\\\
&\texttt{ as.integer(table.old.degrees))}\\\
&\texttt{state\$degree.frequencies \%+\% hset(names(table.new.degrees),}\\\
&\texttt{ as.integer(table.new.degrees))}\\\ \end{split}$ (39)

where state$edge.set and state$degree.frequencies are hset objects containing
the current network and degree distribution, flips$id contains the ties that
are flipped, table.old.degrees and table.new.degrees contain the old and new
degrees for the vertices in $\mathcal{I_{F}}$. Note that the constructor hset
is used to create the set $F$, and the multisets containing the degree
frequencies to be subtracted and added.

In Figure 3 the processes derived from $(X_{t})_{t}$, $({X_{t}}^{-})_{t}$ and
$({X_{t}}^{+})_{t}$ are plotted in red, green and blue respectively. In the
left, the moving average of the acceptance ratio is plotted. In each iteration
the proposed transition is either accepted, or it is not. This binary outcome
is replaced in the plot by the average of 150 binary values around it, giving
an indication on how probable are transitions of state in and out of
equilibrium. The stationary and the sparse chain have a similar behaviour with
transition probability approximately equal to $0.35$ throughout their whole
history. The dense chain starts with a larger transition probability, that is
reduced toward the equilibrium value as the process approaches the stationary
distribution. In the right plot the size of the state, that is the number of
ties in the network is plotted in the three cases, showing how the
distribution of $|{X_{t}}^{-}|$ and $|{X_{t}}^{+}|$ approach the one of
$|{X_{t}}|$ as $t\to\infty$.

Figure 3: Left: acceptance ratio (moving average with $150$ filtered
observations) of the Metropolis-Hastings algorithm. Right: number of ties of
the network (state) at different iterations. Colours distinguish the different
starting points: red for stationary, green for sparse, blue for dense.

In Figure 4 are plotted, for all chains, the empirical cumulative distribution
functions of the degree distribution of states at equally spaced iterations.
For the stationary process $(X_{t},Z_{t})$ these curves, that are plotted in
red, are $\\{\mathrm{ecdf}(Z_{t})\\}$ for $t\in\\{1,1001,2001,...,10001\\}$.
For the processes with sparse and dense starting points, the ECDFs are plotted
in blue and green respectively. The width of the ECDFs is larger when $t$ is
as such, showing how the degree distributions of the networks ${X_{t}}^{-}$
and ${X_{t}}^{+}$ approach the degree distribution of ${X_{t}}$ as
$t\to\infty$.

Figure 4: Empirical cumulative distribution function of the degree
frequencies, at iterations $1$ (lowest line width), $1001$, $2001$, $...$,
$10001$ (largest line width), for the three starting points (stationary - red,
sparse - green, dense - blue).

Markov chains of the type discussed in this application are used to estimate
the parameters of models for which the function that normalizes the stationary
distribution of the process is not analytic. For example in _Exponential
Random Graph Models_ (Robins et al., 2007) computing the normalization
constant requires summing over the set of all possible graphs of a given
dimension, making the computation infeasible also for small networks.
Frequentist and Bayesian estimation algorithms for ERGMs (Snijders et al.,
2002; Caimo and Friel, 2011) are based on the _Markov Chain Monte Carlo
Maximum Liklihood Estimation_ method (_MCMCMLE_), developed in Geyer (1991)
and Geyer and Thompson (1992), where gradients of the likelihood are
approximated using a Markov chain that does not require the computation of the
normalization constant.

A graph is usually represented by its _adjacency matrix_ , with Boolean
elements denoting whether two vertices are connected. Observed large networks
are usually _sparse_ , meaning that the number of edges grows linearly with
respect to the number of vertices. If the adjacency matrix is stored as a
_dense matrix_ , the space required is $O(n^{2})$, where $n$ is the size of
the network, but individual elements can be updated in constant time ($O(1)$
elementary operations). Whereas if a _sparse matrix_ is used, the space
required for storing a sparse network is $O(n)$, but flipping a tie might cost
$O(n)$ elementary operations, as the sparse matrix might have to be re-
constructed. The hash table used in hset can be helpful in these cases, as the
space required to store the set is $O(n)$, but individual elements are updated
with $O(1)$ operations.

##  6 Conclusion

The hset implementation of set operations is motivated by the efficacy of hash
data structure when used with reference semantic, allowing significant
computational advantages in algorithms in which a set or a multiset is used as
container, and it is updated few components at the time. The implementation is
currently restricted to sets or multisets with elements that are numbers (or
sets of numbers), so that mathematical relations between integers and reals
are respected, e.g., that 1L and 1.0 both represents the number $1$. This
approach differs from the library sets where the classes of two objects
determinate whether they can be the same, and reference semantic is not used.

Basic parametrized relations and operations between sets and multisets are
implemented, for most of them (all but the intersection) reference semantic
can speed up algorithms significantly. In R, reference semantic is used for
environment and hash data objects, whereas almost all other objects use value
semantic. Objects with reference semantic are usually modified by functions
with side effects. In Wickham (2019) it is suggested to avoid the use of a
function for both its side effects and its returned value. We partially follow
this suggestion, as for the operations defined in Section 3.2 computed with
reference semantic, the result of the operation is both returned, and the
first operand is transformed to it. Therefore all terms X1%|%X2, X1 <\-
X1%|%X2, and $\texttt{X1 <\- X1\%|}\sim\texttt{\%X2}$ evaluate to the same
objects in memory. The first term is the most common approach to compute
operations with reference semantic, where an operator is viewed as accumulator
causing the first operand to be modified to the result of the operation with
the second operand. However we suggest using the second approach, that is more
coherent with the syntax used with value semantics, and so it is more similar
to how code is usually written in R. Note that the suggestion of avoiding
functions with both side effects and returned value is only partially
followed, because there are no side effects that modify objects that are not
returned.

Our implementation can be useful for simulating Markov chains with countable
state space, as in simulation and estimation of temporal network models.
Recent approaches to statistics and probability theory such as Jacobs (2019),
emphasise the importance of the multiset mathematical structure of discrete
probability distributions, especially when determining the properties of
algorithms used to sample, estimate or learn probability distributions. These
approaches have been heavily influenced by theories of computation, so they
will probably be influential on how algorithms in computational statistics and
other disciplines will be developed and implemented.

## Acknowledgement

Giacomo Ceoldo and Ernst Wit acknowledge funding from the _Swiss National
Science Foundation_ (SNSF 188534).

## References

  * Blitzstein and Diaconis (2011) Joseph Blitzstein and Persi Diaconis.  A sequential importance sampling algorithm for generating random graphs with prescribed degrees.  _Internet mathematics_ , 6(4):489–522, 2011.
  * Brown (2019) Christopher Brown.  _hash: Full Feature Implementation of Hash/Associated Arrays/Dictionaries_ , 2019.  URL https://CRAN.R-project.org/package=hash.  R package version 2.2.6.1.
  * Caimo and Friel (2011) Alberto Caimo and Nial Friel.  Bayesian inference for exponential random graph models.  _Social networks_ , 33(1):41–55, 2011.
  * Geyer (1991) Charles J Geyer.  Markov chain monte carlo maximum likelihood.  _Interface Foundation of North America_ , 1991.
  * Geyer and Thompson (1992) Charles J Geyer and Elizabeth A Thompson.  Constrained monte carlo maximum likelihood for dependent data.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 54(3):657–683, 1992.
  * Holland and Leinhardt (1981) Paul W Holland and Samuel Leinhardt.  An exponential family of probability distributions for directed graphs.  _Journal of the american Statistical association_ , 76(373):33–50, 1981.
  * Jacobs (2019) Bart Jacobs.  _Structured probabilistic reasoning_.  Forthcoming book, 2019.  URL http://www.cs.ru.nl/B.Jacobs/PAPERS/ProbabilisticReasoning.pdf.
  * Meyer and Hornik (2009) David Meyer and Kurt Hornik.  Generalized and customizable sets in R.  _Journal of Statistical Software_ , 31(2):1–27, 2009.  doi: 10.18637/jss.v031.i02.
  * Rinaldo et al. (2013) Alessandro Rinaldo, Sonja Petrović, and Stephen E. Fienberg.  Maximum lilkelihood estimation in the $\beta$-model.  _The Annals of Statistics_ , 41(3):1085 – 1110, 2013.  doi: 10.1214/12-AOS1078.  URL https://doi.org/10.1214/12-AOS1078.
  * Robins et al. (2007) Garry Robins, Pip Pattison, Yuval Kalish, and Dean Lusher.  An introduction to exponential random graph (p*) models for social networks.  _Social networks_ , 29(2):173–191, 2007.
  * Snijders et al. (2002) Tom AB Snijders et al.  Markov chain monte carlo estimation of exponential random graph models.  _Journal of Social Structure_ , 3(2):1–40, 2002\.
  * Wickham (2019) Hadley Wickham.  _Advanced r_.  CRC press, 2019.  URL https://adv-r.hadley.nz/.

##  Appendix A Semantics in R language

In R, objects are generally modified with _value semantic_ , meaning that
whenever the object is accessed, a new copy of the object is first created,
this new copy is then modified, and eventually stored with the name of the
previous variable. More precisely, a new copy is not always created, because
modifying the object locally has computational advantages, this behaviour is
called _copy-on-modified_ and it is explained in Chapter 2 of Wickham (2019),
however when reasoning about the code, it can be assumed that the code behaves
as if an object is copied every time it is accessed. The other approach is to
use a _reference semantic_ , where the name of the object refers to a pointer
to a memory address, so that the modification is always local. However,
operations in the two semantics behave differently. For example, in the code

a = 2; b = a; b = b + 1; (40)

a is initialized to 2, then b is defined to be equal to a, and b is
incremented. With value semantic, after the three operations, a and b are
equal to 2 and 3 respectively, whereas with reference semantic, they are both
equal to 3. The reason is that, with value semantic, in the second operation,
the value that is referred by a is copied and stored elsewhere, this copy is
referred by b, when b is modified, only the value to which b points to is
modified, so a has not changed. On the other hand, if the code above would
have been with reference semantic, the second operation would have copied the
pointer labelled by a, and this new copied pointer is labelled by b, then a
and b would have pointed to the same address in memory, so if b is modified,
also a is.

Objects of class "numeric" (as for most types and classes in R), are accessed
and updated with value semantic. Then after the operations above are
evaluated, a and b are equal to 2 and 3 respectively. Whereas objects of class
"hash" (Brown, 2019), are accessed and updated with reference semantic,
therefore in the code

$\begin{split}&\texttt{a = hash::hash(key="k1", values=2)};\\\ &\texttt{b =
a;}\\\ &\texttt{b[["k1"]] = b[["k1"]] + 1};\end{split}$ (41)

a and b refers to the same hash table, so after the last command, a[["k1"]] is
also equal to 3. An hash table can be copied with the following method

$\begin{split}&\texttt{a = hash::hash(key="k1", values=2)};\\\ &\texttt{b =
hash::hash(hash::keys(a), hash::values(a))};\\\ &\texttt{b[["k1"]] = b[["k1"]]
+ 1};\end{split}$ (42)

so that in the second line a new hash table, labelled by b is created, this
hash table is a copy of a, as the keys and values are copied from it, but the
modification of b in the last row modified, does not change a.