\section{Introduction}

The properties of light vector mesons at finite baryon density, such as their
masses and decay widths, have attracted considerable experimental and
theoretical interest over the last few decades~\cite{Hayano:2008vn,Leupold:2009kz,Hosaka:2016ypm,Krein:2017usp,Metag:2017yuh}. This has been in part
due to their potential to carry information on the partial restoration of
chiral symmetry, and the possible role of QCD of van der Waals forces in the
binding of quarkonia to nuclei. In particular, there is special interest on
the $\phi$-meson, the main reasons being: i) despite its nearly pure
$s\overline{s}$ content, the $\phi$-meson does interact strongly with a
nucleus, composed predominantly of light $u$ and $d$ quarks, through the
excitation of below-threshold virtual kaon and anti-kaon states that might
have their properties changed in medium~\cite{Tsushima:1997df,Laue:1999yv,SchaffnerBielich:1999cp,Akaishi:2002bg,Fuchs:2005zg}; (ii) the $\phi N$ interaction
in vacuum~\cite{Titov:1997qz,Titov:1998bw,Oh:1999nv,Oh:2001bq}and a possible in-medium mass shift of the $\phi$ are related to the strangeness content of the nucleon~\cite{Gubler:2014pta}, which may have
implications beyond the physics of the strong interaction~\cite{Bottino:2001dj,Ellis:2008hf,Giedt:2009mr};
(iii) medium modifications of $\phi$-meson properties have been proposed~\cite{Sibirtsev:2006yk}
as a possible source for the anomalous nuclear mass number $A$-dependence
observed in $\phi$-meson production from nuclear targets~\cite{Ishikawa:2004id};
(iv) furthermore, as the $\phi$-meson is a nearly pure $s\overline{s}$ state
and gluonic interactions are flavor blind, studying it serves to test theories
of the multiple-gluon exchange interactions, including long range QCD van der
Waals forces~\cite{Appelquist:1978rt}, which are believed to play a role in the binding of the
$J/\Psi$ and other exotic heavy-quarkonia to matter~\cite{Brodsky:1989jd,Luke:1992tm,Sibirtsev:1999jr,Gao:2000az,Beane:2014sda,Brambilla:2015rqa,Gao:2017hya,Kawama:2014pja,Kawama:2014iwa,Ohnishi:2014xla,Aoki:2015qla,Morino:2015bqa}. Heavy-ion
collisions and photon- or proton-induced reactions on nuclear targets have
been used to extract information on the in-medium properties of hadrons.
Several experiments have focused on the light vector $\rho$, $\omega$, and
$\phi$ mesons, since their mean-free paths can be comparable with the size
of a nucleus after being produced inside the nucleus. However, a unified
consensus has not yet been reached among the different experiments--see
Refs.~\cite{Hayano:2008vn,Leupold:2009kz,Hosaka:2016ypm,Krein:2017usp,Metag:2017yuh} for comprehensive reviews on the current status.
For the $\phi$-meson, although the precise values are different, a large
in-medium broadening of its decay width has been reported by most of the
experiments performed, while only a few of them find evidence for a
substantial mass shift~\cite{Ishikawa:2004id,Muto:2005za,Mibe:2007aa,Qian:2009ab,Wood:2010ei,Polyanskiy:2010tj}. In 2007 the KEK-E325 collaboration
reported a 3.4\% mass reduction of the $\phi$-meson~\cite{Muto:2005za} and an in-medium
decay width of $\approx$ 14.5 MeV at normal nuclear matter density
$\rho_{0}= 0.15$ fm$^{-3}$. These conclusions were based upon the measurement
of the invariant mass spectra of $e^{+}e^{-}$ pairs in 12 GeV $p+A$ reactions,
with copper and carbon being used as targets~\cite{Muto:2005za}. Even though this result
may indicate a signal for partial restoration of chiral symmetry in nuclear
matter, it is not possible to draw a definite conclusion solely from this.
In fact, recently, a large in-medium $\phi$-meson decay width ($>30$ MeV)
has been extracted at various experimental facilities without observing any
mass shift~\cite{Ishikawa:2004id,Mibe:2007aa,Qian:2009ab,Wood:2010ei,Polyanskiy:2010tj}. It is therefore clear that the search for evidence of
a light vector meson mass shift in nuclear matter is indeed a complicated
issue and further experimental efforts~\cite{Aoki:2015qla,JPARCE16Proposal} are required in order to
understand the phenomenon better. Indeed, the J-PARC E16 collaboration~\cite{Aoki:2015qla,JPARCE16Proposal}
intends to perform a more systematic study for the mass shift of vector mesons
with higher statistics than the above-mentioned experiment at KEK-E325.
However, either complementary or alternative experimental methods are desired.
The study of the $\phi$-meson--nucleus bound states is complementary to the
invariant mass measurements, such as Ref.~\cite{Muto:2005za}, where only a small fraction
of the produced $\phi$-mesons decay inside the nucleus and may be expected to
provide extra information on the $\phi$-meson properties at finite baryon
density. Along these lines, and motivated by the 3.4\% mass reduction reported
by the KEK-E325 experiment~\cite{Muto:2005za}, the E29 collaboration at J-PARC has recently put forward a proposal~\cite{JPARCE29Proposal,JPARCE29ProposalAdd} to study the in-medium mass modification of
the $\phi$-meson via the possible formation of $\phi$-meson--nucleus bound
states~\cite{Ohnishi:2014xla,Buhler:2010zz}. Furthermore, there is also a proposal at JLab, following the
12 GeV upgrade, to study the binding of $\phi$ and $\eta$ mesons to
$^{4}$He~\cite{JLabphiJLabphi}. This new experimental approach~\cite{Ohnishi:2014xla,Aoki:2015qla,Buhler:2010zz,JLabphiJLabphi} for the
measurement of the $\phi$-meson mass shift in nuclei, will produce a slowly
moving $\phi$-meson~\cite{Ohnishi:2014xla,Aoki:2015qla,Buhler:2010zz,JLabphiJLabphi}, where the maximum nuclear matter effect
can be probed. In this way, one may indeed anticipate the formation of a
$\phi$-meson--nucleus bound state, where the $\phi$-meson is trapped inside
the nucleus. Meson-nucleus systems bound by attractive strong interactions
are very interesting objects~\cite{Krein:2017usp,Metag:2017yuh}. First, they are strongly interacting exotic
many-body systems and to study them serves, for example, to understand better
the multiple-gluon exchange interactions, including ``QCD van der Waals''
forces~\cite{Appelquist:1978rt}, which are believed to play a role in the binding of the $J/\Psi$
and other exotic heavy-quarkonia to matter (a nucleus)~\cite{Brodsky:1989jd,Luke:1992tm,Sibirtsev:1999jr,Gao:2000az,Beane:2014sda,Brambilla:2015rqa,Gao:2017hya,Kawama:2014pja,Kawama:2014iwa,Ohnishi:2014xla,Aoki:2015qla,Morino:2015bqa}. Second, they
provide unique laboratories for the study of hadron properties at finite
density, which may not only lead to a deeper understanding of the strong
interaction~\cite{Hayano:2008vn,Leupold:2009kz,Hosaka:2016ypm,Krein:2017usp,Metag:2017yuh} but of the structure of finite nuclei as well~\cite{Stone:2016qmi,Guichon:2006er}.
A downward mass shift of the $\phi$-meson in a nucleus is directly connected
with the possible existence of an attractive potential between the $\phi$-meson
and the nucleus where it has been produced, the strength of which is expected
to be of the same order as that of the mass shift. Concerning the theoretical
evaluation of the $\phi$-meson mass shift, various authors predict a downward
shift of the in-medium $\phi$-meson mass and a broadening of its decay width,
many of them focusing on the self-energy of the $\phi$-meson due to the
kaon-antikaon loop. Ko et al.~\cite{Ko:1992tp} used a density-dependent kaon mass
determined from chiral perturbation theory and found that at normal nuclear
matter density, $\rho_{0}$, the $\phi$-meson mass decreases very little,
by at most 2\%, and a $\Gamma_{\phi}\approx 25$ MeV width which broadens
drastically for large densities. Hatsuda and Lee calculated the in-medium
$\phi$-meson mass based on the QCD sum rule approach~\cite{Hatsuda:1991ez,Hatsuda:1996xt}, and predict a
1.5\%-3\% decrease at $\rho_{0}$. Other investigations also predict a small
downward mass shift and a large broadening of the $\phi$-meson width at
$\rho_{0}$: Ref.~\cite{Klingl:1997tm} reports a negative mass shift of $< 1\%$ and a decay
width of 45 MeV; Ref.~\cite{Oset:2000eg} predicts a decay width of 22 MeV but does not
report a result on the mass shift; and Ref.~\cite{Cabrera:2002hc} gives a rather small
negative mass shift of $\approx 0.81\%$ and a decay width of 30 MeV.
More recently, Ref.~\cite{Gubler:2015yna} reported a downward mass shift of $< 2\%$ and a
large broadening width of 45 MeV at $\rho_{0}$; and finally, in Ref.~\cite{Cabrera:2016rnc},
extending the work of Refs.~\cite{Oset:2000eg,Cabrera:2002hc}, the authors reported a negative mass
shift of 3.4\% and a large decay width of 70 MeV at $\rho_{0}$. The reason
for these differences may lie in the different approaches used to estimate
the kaon-antikaon loop contributions to the $\phi$-meson self-energy and
this might have consequences for the formation of $\phi$-meson--nucleus bound
states. From a practical point of view, the important question is whether
this attraction, if it exists, is sufficient to bind the $\phi$-meson to a
nucleus. A simple argument can be given as follows. One knows that for an
attractive spherical well of radius $R$ and depth $V_{0}$, the condition for
the existence of a nonrelativistic $s$-wave bound state of a particle of mass
$m$ is $V_{0}>\frac{\pi^{2}\hbar^{2}}{8mR^{2}}$. Using $m = m_{\phi}^{*}$, where
$m_{\phi}^{*}$ is the $\phi$-meson mass at normal nuclear matter density found
in Ref.~\cite{Muto:2005za} and $R = 5$ fm (the radius of a heavy nucleus), one obtains
$V_{0}>2$ MeV. Therefore, the prospects of capturing a $\phi$-meson seems
quite favorable, provided that the $\phi$-meson can be produced almost at
rest in the nucleus. Real nuclei, however, have a surface and this estimate
can be quite misleading~\cite{Krein:2017usp}. A full calculation using realistic density
profiles of nuclei is required for a more reliable estimate. The work
presented at this workshop was carried out in collaboration with
professors K.~Tushima, G.~Krein, and A.~W.~Thomas and has been published in
Refs.~\cite{Cobos-Martinez:2017vtr,Cobos-Martinez:2017woo}. In Ref.~\cite{Cobos-Martinez:2017vtr} we studied the $\phi$-meson mass shift and
decay width in nuclear matter, based on an effective Lagrangian approach, by
evaluating the $K\overline{K}$ loop contribution in the $\phi$-meson
self-energy, with the in-medium K and $\overline{K}$ masses explicitly
calculated by the quark-meson coupling (QMC) model~\cite{Saito:2005rv}. This initial study
has been extended in Ref.~\cite{Cobos-Martinez:2017woo} to some selected nuclei by computing the
$\phi$-meson--nucleus bound complex potential in the local density
approximation and solving the Klein-Gordon equation in order to obtain the
bound state energies and absorption widths. The nuclear density distributions
for all nuclei studied, except for $^{4}$He, are explicitly calculated using the
QMC model~\cite{Saito:1996sf}.

\section{$\phi$-meson mass and decay width in nuclear matter}

The $\phi$-meson property modifications in nuclear matter, such as its mass
and decay width, are strongly correlated to its coupling to the
$K\overline{K}$ channel, which is the dominant decay channel in vacuum.
Thus, one expects that a significant fraction of the density dependence of
the $\phi$-meson self-energy in nuclear matter arises from the in-medium
modification of the KK intermediate state in the $\phi$-meson self-energy.
We briefly review the computation~\cite{Cobos-Martinez:2017vtr} of the $\phi$-meson self-energy in
vacuum and in nuclear matter using an effective Lagrangian approach~\cite{Klingl:1996by}.
The interaction Lagrangian $\mathcal{L}_{\text{int}}$ involves
$\phi K\overline{K}$ and $\phi\phi K\overline{K}$ couplings dictated by a
local gauge symmetry principle:
\begin{equation}
  \label{eqn:Lint}
  \mathcal{L}_{\text{int}}=  \mathcal{L}_{\phi K\overline{K}}
  +\mathcal{L}_{\phi\phi K\overline{K}}
  \end{equation}
\noindent where (we use the convention
$K=\left(\begin{array}{c} K^{+} \\ K^{0} \end{array}\right)$ for isospin
doublets)
\begin{eqnarray}
\label{eqn:Lpkk}
\mathcal{L}_{\phi K\overline{K}}&=&\mi g_{\phi}\phi^{\mu}
\left[\overline{K}(\partial_{\mu}K)-(\partial_{\mu}\overline{K})K\right], \\
\label{eqn:Lppkk}
\mathcal{L}_{\phi\phi K\overline{K}}&=& g^2_{\phi} \phi^\mu\phi_\mu \overline{K} K.
\end{eqnarray}
We note that the use of the effective interaction Lagrangian of \eqn{eqn:Lint}
without the term given in \eqn{eqn:Lppkk} may be considered as being motivated
by the hidden gauge approach in which there are no four-point vertices,
such as \eqn{eqn:Lppkk}, that involve two pseudoscalar mesons and two vector
mesons~\cite{Lin:1999ve,Lee:1994wx}. This is in contrast to the approach of using the minimal
substitution to introduce vector mesons as gauge particles where such
four-point vertices do appear. However, these two methods have been shown to
be consistent if both the vector and axial vector mesons are included~\cite{Yamawaki:1986zz,Meissner:1986tc,Meissner:1987ge,Saito:1987ba}.
Therefore, we present results with and without such an interaction~\cite{Cobos-Martinez:2017vtr,Cobos-Martinez:2017woo}.
We consider first the contribution from the $\phi K\overline{K}$ coupling,
given by \eqn{eqn:Lpkk}, to the scalar part of the $\phi$-meson self-energy,
$\Pi_{\phi}(p)$.
\begin{equation}
\label{eqn:phise}
\mi\Pi_{\phi}(p)=-\frac{8}{3}g_{\phi}^{2}\int\dfd{4}{4}{q}\vec{q}^{\,2}
D_{K}(q)D_{K}(q-p) \, ,
\end{equation}
\noindent where $D_{K}(q)=\left(q^{2}-m_{K}^{2}+\mi\epsilon\right)^{-1}$ is the
kaon propagator;  $p=(p^{0}=m_{\phi},\vec{0})$ is the $\phi$ meson four-momentum
vector at rest, with $m_{\phi}$ the $\phi$ meson mass; $m_{K} (=m_{\overline{K}})$ is
the kaon mass, and $g_{\phi}$ is the coupling constant. When $m_{\phi}<2m_{K}$
the self-energy $\Pi_{\phi}(p)$ is real. However, when  $m_{\phi}>2m_{K}$, which
is the case here, $\Pi_{\phi}(p)$ acquires an imaginary part. The integral in
\eqn{eqn:phise} is divergent but it will be regulated using a phenomenological
form factor, with cutoff parameter $\Lambda_{K}$, as in Refs.~\cite{Cobos-Martinez:2017vtr,Krein:2010vp}.
The sensitivity of the results to the cutoff value is analyzed below. The mass and
decay width of the $\phi$-meson in vacuum ($m_{\phi}$ and $\Gamma_{\phi}$), as
well as in nuclear matter ($m_{\phi}^{*}$ and $\Gamma_{\phi}^{*}$), are
determined self-consistently~\cite{Cobos-Martinez:2017vtr} from
\begin{eqnarray}
\label{eqn:phimass}
m_{\phi}^{2}&=&\left(m_{\phi}^{0}\right)^{2}+\Re\Pi_{\phi}(m_{\phi}^{2}), \\
\label{eqn:phiwidth}
\Gamma_{\phi}&=&-\frac{1}{m_{\phi}}\Im\Pi_{\phi}(m_{\phi}^{2}).
\end{eqnarray}
The coupling constant $g_{\phi}$ is determined~\cite{Cobos-Martinez:2017vtr} from \eqns{eqn:phimass}
and~(\ref{eqn:phiwidth}), and the experimental value for the $\phi\to K\overline{K}$
decay width in vacuum, corresponding to the branching ratio of 83.1\% of the
total decay width~\cite{PDG:2015}.
The nuclear density dependence of the $\phi$-meson mass and decay width is
driven by the intermediate $K\overline{K}$ state interactions with the nuclear medium.
This effect enters through $m_{K}^{*}$ in the kaon propagators in
\eqn{eqn:phise}. The in-medium mass $m_{K}^{*}$is calculated within the QMC
model~\cite{Cobos-Martinez:2017vtr}, which has proven to be very successful in studying the properties
of hadrons in nuclear matter and finite nuclei. For a more complete discussion
of the model see Refs.~\cite{Tsushima:1997df,Saito:2005rv,Guichon:1989tx}. Here we just make a few necessary comments.
In order to calculate the in-medium properties of $K$ and $\overline{K}$, we consider
infinitely large, uniformly symmetric, spin-isospin-saturated nuclear matter
in its rest frame, where all the scalar and vector mean field potentials,
which are responsible for the nuclear many-body interactions, become constant
in the Hartree approximation~\cite{Cobos-Martinez:2017vtr}.
In \fig{fig:nuclmatt} (left panel) we present the resulting in-medium kaon
Lorentz scalar mass (=antikaon Lorentz scalar mass), calculated using the QMC
model, as a function of the baryon density. The kaon effective mass at normal
nuclear matter density $\rho_{0}= 0.15$ fm$^{-3}$ has decreased by about 13\%.
We also recall, in connection with the calculation of the in-medium $K\overline{K}$
loop contributions to the $\phi$-meson self-energy, that the isoscalar-vector
$\omega$ mean field potentials arise both for the kaon and antikaon. However,
they have opposite signs and cancel each other. Equivalently, they can be
eliminated by a variable shift in the loop calculation~\cite{Tsushima:1997df,Saito:2005rv,Guichon:1989tx} of the
$\phi$-meson self-energy, and therefore we do not show them here.
\begin{figure}
  \centering
\scalebox{0.9}{
\begin{tabular}{ccc}
 \includegraphics[scale=0.2]{mK_in_medium} &
 \includegraphics[scale=0.2]{Dmphi_in_medium} &
 \includegraphics[scale=0.176]{Gphi_in_medium} \\
\end{tabular}
}
\caption{\label{fig:nuclmatt} Left panel: In-medium kaon mass (=antikaon)
  Lorentz scalar mass $m_{K}^{*}$; centre and right panels: $\phi$-meson mass
  shift and decay width, respectively, in symmetric nuclear matter for three
  values of the cutoff parameter $\Lambda_{K}$.}
\end{figure}
In \fig{fig:nuclmatt}, we present the $\phi$-meson mass shift (centre panel)
and decay width (right panel) as a function of the nuclear matter density,
$\rho_{B}$, for three values of the cutoff parameter $\Lambda_{K}$. As can be
seen, the effect of the in-medium kaon and antikaon mass change yields a
negative mass shift for the $\phi$-meson. This is because the reduction in the
kaon and antikaon masses enhances the $K\overline{K}$ loop contribution in nuclear
matter relative to that in vacuum. For the largest value of the nuclear matter
density, the downward mass shift turns out to be a few percent at most for
all values of $\Lambda_{K}$. On the other hand, we see that $\Gamma_{\phi}^{*}$
is very sensitive to the change in the kaon and antikaon masses, increasing
rapidly with increasing nuclear matter density, up to a factor of $\approx 20$
enhancement for the largest value of $\rho_{B}$. In \tab{tab:nuclmatt} we
present the values for $m_{\phi}^{*}$and $\Gamma_{\phi}^{*}$ at normal nuclear
matter density $\rho_{0}$. We see that the negative kaon and antikaon mass
shift of 13\%~\cite{Cobos-Martinez:2017vtr} induces a downward mass shift of the $\phi$-meson of just
$\approx 2\%$, while the broadening of the $\phi$-meson decay width is an
order-of-magnitude larger than its vacuum value.
\begin{table}
  \centering
\scalebox{1.0}{
  \begin{tabular}{c|ccc}
    \hline
    & $\Lambda_{K}=2000$ & $\Lambda_{K}= 3000$ & $\Lambda_{K}= 4000$ \\
    \hline
$m_{\phi}^{*}$ & 1000.9 (1009.5) & 994.9 (1004.3) & 990.4 (1000.6) \\
$\Gamma_{\phi}^{*}$ & 34.8 (37.8) & 32.8 (36.0) & 31.3 (34.7) \\
  \hline
\end{tabular}
}
\caption{\label{tab:nuclmatt} $\phi$-meson mass and width at normal nuclear
  matter density $\rho_{0}$ with and without the gauged Lagrangian of
  \eqn{eqn:Lppkk}. The values in parentheses were computed by adding the
  gauged Lagrangian of \eqn{eqn:Lppkk}. All quantities are given in MeV.}
\end{table}
In Ref.~\cite{Cobos-Martinez:2017vtr} we evaluated the impact of adding the $\phi\phi K\overline{K}$
interaction of \eqn{eqn:Lppkk} on the in-medium $\phi$-meson mass and
decay width. We found that one still gets a downward shift of the in-medium
$\phi$-meson mass as well as a significant broadening of the decay width when
this interaction is added. In both cases, though, the absolute values are
slightly different from those shown in \fig{fig:nuclmatt}.
In tab{tab:nuclmatt} we present the values for $m_{\phi}^{*}$ and
$\Gamma_{\phi}^{*}$ at $\rho_{0}$ obtained by adding the gauged Lagrangian of
\eqn{eqn:Lppkk}. In both cases, for the mass and decay width in nuclear
matter, the effect of adding \eqn{eqn:Lppkk} can be compensated by the use
of a larger cutoff $\Lambda_{K}$.
The results described above support those which suggest that one should
observe a small downward mass shift and a large broadening of the decay width
of the $\phi$-meson in a nuclear medium. Furthermore, they open experimental
possibilities for studying the binding and absorption of $\phi$-meson in
nuclei. Although the mass shift found in this study may be large enough to
bind the $\phi$-meson to a nucleus, the broadening of its decay width will
make it difficult to observe a signal for the $\phi$-meson--nucleus bound
state formation experimentally. We explore this further in the second part of
this talk.

\section{$\phi$-meson--nucleus bound states }

In this part we discuss the situation where the $\phi$-meson is placed in a
nucleus. The nuclear density distributions for the nuclei $^{12}$C, $^{16}$O,
and $^{208}$Pb are obtained using the QMC model~\cite{Saito:1996sf}. For $^{4}$He, we use the
parametrization for the density distribution obtained in Ref.~\cite{Saito:1997ae}. Then,
using a local density approximation we calculate the $\phi$-meson complex
potentials for a nucleus $A$, which can be written as
\begin{equation}
\label{eqn:Vcomplex}
V_{\phi A}(r)= U_{\phi}(r)-\frac{\mi}{2}W_{\phi}(r),
\end{equation}
\noindent where $r$ is the distance from the center of the nucleus and
$U_{\phi}(r)=\Delta m_{\phi}(\rho_{B}(r)) \equiv m^{*}_{\phi}(\rho_{B}(r))-m_\phi$
and $W_{\phi}(r)=\Gamma_{\phi}(\rho_{B}(r))$ are, respectively,
the $\phi$-meson mass shift and decay width in a nucleus $A$, with
$\rho_{B}(r)$ the baryon density distribution for the particular nucleus.
\begin{figure}
  \centering
\scalebox{0.675}{
\begin{tabular}{cccc}
  \includegraphics[scale=0.2]{Vphi_he4} &
  \includegraphics[scale=0.2]{Gphi_he4} &
  \includegraphics[scale=0.2]{Vphi_c12} &
  \includegraphics[scale=0.2]{Gphi_c12} \\
  \includegraphics[scale=0.2]{Vphi_o16} &
  \includegraphics[scale=0.2]{Gphi_o16} &
  \includegraphics[scale=0.2]{Vphi_pb208} &
  \includegraphics[scale=0.2]{Gphi_pb208} \\
\end{tabular}
}
\caption{\label{fig:phinuclpot} Real [$U_{\phi}(r)$] and imaginary [$W\phi(r)$]
  parts of the $\phi$-meson--nucleus potentials for the nuclei $^{4}$He,
  $^{12}$C, $^{16}$O, and $^{208}$Pb, for three values of the cutoff parameter
  $\Lambda_{K}$.}
\end{figure}
\fig{fig:phinuclpot} shows the $\phi$-meson potentials calculated for the
selected nuclei, for three values of the cutoff parameter $\Lambda_{K}$. One
can see that the depth of the real part of the potential is sensitive to the
cutoff parameter, varying from -20 MeV to -35 MeV for $^{4}$He and from
-20 MeV to -30 MeV for $^{208}$Pb. On the other hand, one can see that the
imaginary part does not vary much with $\Lambda_{K}$. These latter observations
may well have consequences for the feasibility of experimental observation
of the expected bound states.
Using the $\phi$-meson potentials obtained in this manner, we next calculate
the $\phi$-meson--nucleus bound state energies and absorption widths for the
selected nuclei. Before proceeding, a few comments on the use of \eqn{eqn:kge}
are in order. In this study, we consider the situation where the $\phi$-meson
is produced nearly at rest. Then, it should be a very good approximation to
neglect the possible energy difference between the longitudinal and transverse
components of the $\phi$-meson wave function $\psi_{\phi}^{\mu}$. After imposing
the Lorentz condition, $\partial_{\mu}\psi_{\phi}^{\mu}=0$, to solve the Proca
equation becomes equivalent to solving the Klein-Gordon equation
\begin{equation}
\label{eqn:kge}
\left(-\nabla^{2} + \mu^{2} + 2\mu V(\vec{r})\right)\phi(\vec{r})
= \mathcal{E}^{2}\phi(\vec{r}),
\end{equation}
where $μ = m_{\phi}m_{A}/(m_{\phi} + m_{A})$ is the reduced mass of the
$\phi$-meson--nucleus system with $m_{\phi}$ ($m_{A}$) the mass of the
$\phi$-meson (nucleus $A$) in vacuum, and $V (r)$ is the complex
$\phi$-meson--nucleus potential of \eqn{eqn:Vcomplex}. We solve the
Klein-Gordon equation using the momentum space methods~\cite{Kwan:1978zh}.
The calculated bound state energies ($E$) and absorption widths ($\Gamma$), which are
related to the complex energy eigenvalue $\mathcal{E}$ by
$E= \Re\mathcal{E}-\mu$ and $\Gamma= -2\Im\mathcal{E}$, are listed in
\tab{tab:phienergies} for three values of the cutoff parameter $\Lambda_{K}$,
with and without the imaginary part of the potential.
We first discuss the case in which the imaginary part of the
$\phi$-meson--nucleus potential is set to zero. The results are given within
parentheses in \tab{tab:phienergies}. From the values shown, we see that the
$\phi$-meson is expected to form bound states with all the nuclei selected,
for all values of the cutoff parameter $\Lambda_{K}$. The bound state energy
is obviously dependent on $\Lambda_{K}$, increasing as $\Lambda_{K}$ increases.
Next, we discuss the results obtained when the imaginary part of the potential
is included. Adding the absorptive part of the potential changes the situation
appreciably. From the results presented in \tab{tab:phienergies} we note that
for the largest value of the cutoff parameter $\Lambda_{K}$, which yields the
deepest attractive potentials (see \fig{fig:phinuclpot}), the $\phi$-meson
is expected to form bound states with all the selected nuclei, including the
lightest one, the $^{4}$He nucleus. However, in this case, whether or not the
bound states can be observed experimentally is sensitive to the value of the
cutoff parameter $\Lambda_{K}$. One also observes that the width of the bound
state is insensitive to the values of $\Lambda_{K}$ for all nuclei.
Furthermore, since the so-called dispersive effect of the absorptive potential
is repulsive, the bound states disappear completely in some cases, even though
they were found when the absorptive part was set to zero. This feature is
obvious for the $^{4}$He nucleus, making it especially relevant to the future
experiments, planned at J-PARC and JLab using light and medium-heavy
nuclei~\cite{Ohnishi:2014xla,Aoki:2015qla,Buhler:2010zz,JLabphiJLabphi}.
\begin{table}[ht]
\begin{center}
\scalebox{0.85}{
  \begin{tabular}{ll|rr|rr|rr}
\hline \hline
& & \multicolumn{2}{c|}{$\Lambda_{K}=2000$} &
\multicolumn{2}{c}{$\Lambda_{K}=3000$} &
\multicolumn{2}{c}{$\Lambda_{K}=4000$}  \\
\hline
 & & $E$ & $\Gamma/2$ & $E$ & $\Gamma/2$ & $E$ & $\Gamma/2$ \\
\hline
$^{4}_{\phi}\text{He}$ & 1s & n (-0.8) & n & n (-1.4) & n & -1.0 (-3.2) & 8.3 \\
\hline
$^{12}_{\phi}\text{C}$ & 1s & -2.1 (-4.2) & 10.6 & -6.4 (-7.7) & 11.1 & -9.8
(-10.7) & 11.2 \\
\hline
$^{16}_{\phi}\text{O}$ & 1s & -4.0 (-5.9) & 12.3 & -8.9 (-10.0) & 12.5 & -12.6
(-13.4) & 12.4 \\
& 1p & n (n) & n & n (n) & n & n (-1.5) & n \\
\hline
$^{208}_{\phi}\text{Pb}$ & 1s & -15.0 (-15.5) & 17.4 & -21.1 (-21.4) & 16.6 &
-25.8 (-26.0) & 16.0 \\
& 1p & -11.4 (-12.1) & 16.7 & -17.4 (-17.8) & 16.0 & -21.9  (-22.2) & 15.5 \\
& 1d & -6.9 (-8.1) & 15.7 & -12.7 (-13.4) & 15.2 & -17.1 (-17.6) & 14.8 \\
& 2s & -5.2 (-6.6) & 15.1 & -10.9 (-11.7) & 14.8 & -15.2 (-15.8) & 14.5 \\
& 2p & n (-1.9) & n & -4.8 (-6.1) & 13.5 & -8.9 (-9.8) & 13.4 \\
& 2d & n (n) & n & n (-0.7) & n & -2.2 (-3.7) & 11.9 \\
\hline \hline
\end{tabular}}
\caption{\label{tab:phienergies} $\phi$-nucleus single-particle energies $E$
  and half width $\Gamma/2$ obtained, with and without the imaginary part of
  the potential of \eqn{eqn:Vcomplex}, for three values of the cutoff parameter
  $\Lambda_K$. When only the real part is included, the corresponding
  single-particle energy $E$ is given inside parenthesis and $\Gamma=0$ for
  all nuclei. ``n'' indicates that no bound state is found. All quantities are
  given in MeV.}
\end{center}
\end{table}

\section{Summary and discussion}

We have calculated the $\phi$-meson mass and width in nuclear matter within an
effective Lagrangian approach up to three times of normal nuclear matter
density. Essential to our results are the in-medium kaon masses, which are
calculated in the QMC model, where the scalar and vector meson mean fields
couple directly to the light $u$ and $d$ quarks (antiquarks) in the $K$
($\overline{K}$) mesons. At normal nuclear matter density, allowing for a very large
variation of the cutoff parameter $\Lambda_{K}$, although we have found a
sizable negative mass shift of 13\% in the kaon mass, this induces only a
few percent downward shift of the $\phi$-meson mass. On the other hand, it
induces an order-of-magnitude broadening of the decay width. We have also
calculated the $\phi$-meson–nucleus bound state energies and absorption widths
for various nuclei. The $\phi$-meson–nucleus potentials were calculated using a
local density approximation, with the inclusion of the $K\overline{K}$ loop in the
$\phi$-meson self-energy. The nuclear density distributions, as well as the
in-medium $K$ and $\overline{K}$ meson masses, were consistently calculated by
employing the quark-meson coupling model. Using the $\phi$-meson--nucleus
complex potentials found, we have solved the Klein-Gordon equation in momentum
space, and obtained $\phi$-meson--nucleus bound state energies and absorption
widths. Furthermore, we have studied the sensitivity of our results to the
cutoff parameter $\Lambda_{K}$ in the form factor at the $\phi K\overline{K}$ vertex
appearing in the $\phi$-meson self-energy. We expect that the $\phi$-meson
should form bound states for all four nuclei selected, provided that the
$\phi$-meson is produced in (nearly) recoilless kinematics. This feature is
even more obvious in the (artificial) case where the absorptive part of the
potential is ignored. Given the similarity of the binding energies and widths
reported here, the signal for the formation of the $\phi$-meson--nucleus bound
states may be difficult to identify experimentally. Therefore, the feasibility
of observation of the $\phi$-meson–nucleus bound states needs further
investigation, including explicit reaction cross section estimates.

\section{Acknowledgements}

This work was partially supported by Conselho Nacional de Desenvolvimento
Cient\'{\i}fico e Tecnol\'ogico-CNPq, Grants No.~152348/2016-6 (J.J.C-M.),
No.~400826/2014-3 and No.~308088/2015-8 (K.T.), No.~305894/2009-9 (G.K.), and
No.~313800/2014-6 (A.W.T.), and Funda{\c c}\~{a}o de Amparo \`{a} Pesquisa do
Estado de S\~ao Paulo-FAPESP, Grants No.~2015/17234-0 (K.T.) and
No.~2013/01907-0 (G.K.). This research was also supported by the University of
Adelaide and by the Australian Research Council through the ARC Centre of
Excellence for Particle Physics at the Terascale (CE110001104), and through
Grant No.~DP151103101 (A.W.T.). J.J.C-M. also acknowledges the support given
by Red-FAE-CONACyT (M\'exico) to attend the workshop where this work was
presented.

\section*{References}


\section{Signal Temporal Logic}
\label{sec:background}

\paragraph*{ \bf STL} {\it Signal Temporal Logic} (STL)~\cite{Maler2004} is a formal
specification language to express temporal properties over
real-values trajectories with dense-time interval. For the rest of the paper, let be $\vec x: \mathbb{T} \rightarrow \mathbb{R}^n$ a trace/trajectory, where $ \mathbb{T}  = \mathbb{R}_{\geq 0}$ is the time domain, $x_i(t)$ is the value at time $t$ of the projection on the $i^{th}$ coordinate, and $\vec x = (x_1, ..., x_n)$, as an abuse on the notation, is used also to indicate the set of variables of the trace considered in the formulae.
\begin{definition} [STL syntax] The syntax of STL is given by
$$\varphi :=
    \top \:|\: \mu  \:|\:  \neg \varphi \:|\: \varphi_{1} \wedge \varphi_{2}  \:|\:  \varphi_{1} \: \until{I} \: \varphi_{2}  $$
where  $\top$ is the Boolean true constant, $\mu$ is an {\it atomic proposition}, inequality of the form $y(x) > 0$ ($y:\mathbb{R}^n \rightarrow \mathbb{R}$), {\it negation} $\neg$ and {\it conjunction} $\wedge$ are the standard Boolean connectives, and $\until{I}$ is the {\it Until}  temporal modality, where  $I$ is a real positive interval.
As customary, we can derive the {\it disjunction} operator $\vee$ and the future  {\it eventually} $\ev{I}$ and  {\it always} $\glob{I}$ operators from  the until temporal modality $( \varphi_1 \vee \varphi_2 = \neg (\neg \varphi_1 \wedge \neg \varphi_2),$ $\ev{I} \varphi = \top \until{I} \varphi,\mbox{ } \glob{I} \varphi = \neg  \ev{I} \neg \varphi ).$
\end{definition}

 STL can be interpreted over a trajectory $\vec x$ using a
 qualitative (yes/no) or a quantitative (real-value) semantics~\cite{Maler2004,Donze2013}.

 We report here only the quantitative semantics and we refer the reader
 to~\cite{MalerN13,Maler2004,Donze2013} for more details.
\begin{definition} [{STL  Quantitative Semantics}]
\label{sp.rob}
The \emph{quantitative satisfaction function} $\rho$ returns a value
$\rho(\varphi, \vec{x},t) \in \bar{\mathbb{R}}$,\footnote{$\bar{\mathbb{R}}= \mathbb{R}\cup \{ -\infty, +\infty\}.$}
quantifying the robustness degree (or satisfaction degree) of the property $\varphi$ by the trajectory $\vec{x}$
at time $t$.  It is defined recursively as follows:
{\small
\begin{align*}
&\rho(\top, \vec x, t)  & = &\mbox{ } \phantom{a}  {\large + \infty} \\
&\rho(\mu, \vec{x},t)  & = &\mbox{ } \phantom{a}  y(\vec{x}(t))  \mbox{   where } \mu \equiv y(\vec{x}(t)) \geq 0 \\
& \rho (\neg \varphi,\vec{x},t)  & =  & \mbox{ } \phantom{a}  - \rho (\varphi,\vec{x},t)\\
&\rho( \varphi_{1} \wedge  \varphi_{2}, \vec{x},t)  & = &\mbox{ } \phantom{a} \min ( \rho( \varphi_{1},\vec{x},t),\rho( \varphi_{2},\vec{x},t) ) \\
& \rho(  \varphi_{1} \: \until{[a,b]}  \varphi_{2}, \vec{x},t) & = &\mbox{ } \phantom{a}  \sup_{t'\in t+[a,b]}(\min(\rho( \varphi_{2},\vec{x},t'),\inf_{t'' \in [t,t')}(\rho( \varphi_{1},\vec{x},t''))))
\end{align*}
}
Moreover, we let $\rho(\varphi, \vec{x}):=\rho(\varphi, \vec{x},0)$.
\end{definition}
The sign of $\rho(\varphi, \vec{x})$ provides the link with the standard Boolean semantics
of~\cite{Maler2004}: $\rho(\varphi, \vec{x})>0$  only if  $\vec{x}\models \varphi$, while $\rho(\varphi, \vec{x})<0$
only if $\vec{x}\not\models\varphi$\footnote{The case  $\rho(\varphi, \vec{x}) = 0$, instead, is a borderline case, and
the truth of $\varphi$ cannot be assessed from the robustness degree alone.}.
The absolute value of $\rho(\varphi, \vec{x})$, instead,  can be interpreted as a measure of the robustness of the satisfaction with respect to noise in signal $\vec{x}$, measured in terms of the maximal
perturbation in the secondary signal $  y(\vec{x}(t)) $, preserving truth value. This means that if $\rho(\varphi,\vec{x},t)=r$  then for every signal
$\vec{x}'$ for which every secondary signal satisfies $\max_{t}| y_{j}(t) - y'_{j}(t)|<r$, we have that
$\vec{x}(t) \models \varphi$ if and only if $\vec{x'}(t)\models\varphi$ ({\it correctness property}).

\paragraph*{\bf PSTL} {\it Parametric Signal Temporal Logic}~\cite{AsarinDMN11} is an extension
of STL that parametrizes the formulas. There are two types of formula parameters:
temporal parameters, corresponding to the time bounds in the time intervals associated
with temporal operators (e.g. $a,b \in \mathbb{R}_{\geq 0}$, with $a<b$, s.t. $\ev{[a,b]} \mu$),
and the threshold parameters, corresponding to the constants in the inequality predicates
(e.g., $k \in \mathbb{R}, \mu = x_i > k$, where $x_i$ is a variable of the trajectory).
In this paper, we allow only atomic propositions of the form $\mu = x_i \bowtie k$ with
 $\bowtie \in \{ >, \leq \}$.

Given an STL formula $\phi$, let  $ \mathbb{K}=(\mathcal{T}\times\mathcal{C})$ be
the {\it parameter space}, where $\mathcal{T} \subseteq \mathbb{R}_{\geq 0}^{n_t}$
is the temporal parameter space ($n_t$ being the number of temporal parameters),
and $\mathcal{C}   \subseteq  \mathbb{R}^{n_k}$ is the threshold parameter space
($n_k$ being the number of threshold parameters). A $\boldsymbol{\theta} \in \mathbb{K}$
is a parameter configuration that induces a corresponding STL formula; e.g.,
$\phi=\ev{[a,b]}( x_i > k), \boldsymbol{\theta} = (0,2,3.5)$ then $\phi_{ \boldsymbol{\theta}}=\ev{[0,2]}( x_i > 3.5)$.

\paragraph*{\bf  Stochastic Robustness} Let us consider an unknown stochastic process $({\bf X}(t))_{t \in T}=(X_{1}(t),...,X_{n}(t))_{t \in T}$, where each vector ${\vec X}(t)$ corresponds to the state of the system at time $t$.
For simplicity, we indicate the stochastic model with $\vec X(t)$.
$\vec{X}(t)$ is a measurable also as a random variable $\vec{X}$ on the space $\sspace{}$-valued {\it cadlag functions} ${\cal{D}}([0,\infty),\sspace{})$, here denoted by ${\cal{D}}$, assuming the domain $\pstsp$ to be fixed. It means that the set of  trajectories $\vec x$  of the stochastic process $\vec{X}$ is the set ${\cal{D}}$. Let us consider now an STL formula $\phi$, with predicates interpreted over state variables of ${\vec X}$. Given a trajectory $\vec{x}(t)$ of a stochastic system,  its robustness $\rho(\phi,\vec{x},0)$ is a measurable functional $R_\phi$~\cite{BartocciBNS15} from the trajectories in $\cal{D}$ to $\mathbb{R}$
which defines the real-valued random variable $R_{\phi} = R_{\phi}(\vec{X})$ with probability distribution:
 \[\mathbb{P}\left(R_{\phi}(\vec{X})\in [a,b] \right)= \mathbb{P}\left( \vec{X} \in \{\vec{x}\in\mathcal{D}~|~\rho(\phi,\vec{x},0)\in [a,b] \}   \right).  \]
Such distribution of robustness degrees can be summarized by the average robustness degree. Fixing the stochastic process $\vec X$, $\mathbb{E}(R_{\phi}| \vec{X})$,  it gives a measure of how strongly a formula is satisfied.
The satisfaction is more robust when this value is higher. In this paper, we will approximate this expectation by Monte Carlo sampling.

\section{Conclusion}
\label{sec:conclusion}

We present a framework to learn from a labeled dataset of normal and anomalous trajectories a {\it Signal Temporal Logic} (STL)
specification
that better discriminates among them.
In particular, we design a Robust Genetic algorithm (ROGE) that
combines an {\it Evolutionary Algorithm} (EA) for learning
the structure of the  formula  and the {\it Gaussian Process Upper Confidence Bound}
algorithm (GP-UCB) for synthesizing the formula parameters.
We compare ROGE with our previous work \cite{BufoBSBLB14}
and with the Decision Tree approach presented in~\cite{bombara_decision_2016} on
an anomalous trajectory detection problem of a maritime surveillance system and on an Ineffective Inspiratory Effort example.

A significant difference concerning our previous approach~\cite{BBS14, BufoBSBLB14} is that we avoid  reconstructing a generative statistical model from the dataset. Furthermore, we modified both the structure and parameter optimization procedure of the genetic algorithm, which now relays on the robustness semantics of STL. This structural improvement was necessary considering that a naive application of our previous approach~\cite{BufoBSBLB14}
directly on datasets performs very poorly.

We compare our method also with the Decision Tree (DTL4STL) approach
of~\cite{bombara_decision_2016} obtaining a low misclassification rate and producing smaller and more interpretable STL specifications. Furthermore,  we do not restrict the class of the temporal formula to only \emph{eventually} and \emph{globally} and we allow arbitrary temporal nesting.

Our genetic algorithm can get wrong results if the formulae of the initial generation are chosen entirely randomly. We avoid this behavior by considering simpler formulae from the beginning as a result of the \textsc{generateInitialFormulae} routine. This approach is a form of regularization and resembles the choice of the set of primitive of DTL4STL.

As future work, we plan to develop an iterative method which uses the proposed genetic algorithm and the robustness of STL to reduce the misclassification rate of the generated formula.
The idea is to use the robustness value of a learned formula $\phi$ to identify the region of the trajectory space $\mathcal{D}$ where the generated formula $\phi$ has a high accuracy, i.e. trajectories whose robustness is greater than a positive threshold $h^+$ or smaller than a negative threshold $h^-$. These thresholds can be identified by reducing the number of false positives or false negatives. We can then train a new STL classifier $\phi'$ on the remaining trajectories, having small robustness for $\phi$. This will create a hierarchical classifier, that first tests on $\phi$,  if robustness is too low it tests on $\phi'$, and so on. The depth of such classification is limited only by the remaining data at each step. The method can be further strengthen by relying on bagging to generate an ensemble of classifiers at each level, averaging their predictions or choosing the best answer, i.e. the one with larger robustness.

\section{Introduction}
\label{sec:introduction}
Learning temporal logic requirements from data is an emergent
research field gaining momentum in the rigorous engineering of
cyber-physical systems.
Classical machine learning methods typically generate
very powerful black-box (statistical) models. However,
these models often do not help in the comprehension of the
phenomenon they capture.
Temporal logic provides a precise formal specification language
that can easily be  interpreted by humans.
The possibility to describe datasets in a concise way using temporal logic formulas can thus help to better clarify and comprehend which are the emergent patterns for
the system at hand.
A clearcut example is the problem of anomaly
detection, where the input is a set of trajectories describing
regular or anomalous behaviors, and the goal is to learn a
classifier that not only can be used to detect anomalous
behaviors at runtime, but also gives insights on what
characterizes an anomalous behavior.
Learning temporal properties is also relevant in combination
with state of the art techniques for search-based falsification
of complex closed-loop systems~\cite{0001SDKJ15,Sankaranarayanan17,BartocciBNS15,SilvettiPB17},
as it can provide an automatic way to describe desired
(or  unwanted behaviors) that the system needs to satisfy.

In this paper, we consider the problem of learning a temporal
logic specification from a set of trajectories which
are labeled by human experts (or by any other method
which is not usable for real-time monitoring)
as ``good" for the normal expected behaviors and
``bad" for the anomalous ones.
The goal is to automatically synthesize both
the structure of the formula and its parameters
providing a temporal logic classifier that
can discriminate as much as possible the
bad and the good behaviors.
This specification can be turned into
a monitor that output a positive verdict
for good behaviors and a negative verdict
for bad ones.

\paragraph*{\textbf{Related Work}}
\input{relatedwork}

\paragraph*{\textbf{Our contribution}}

In this work, we consider the problem of mining the formula directly from a data set without requiring to learn a generative model from data.
To achieve this goal,  we introduce a number of techniques to improve the potentials of the genetic algorithm and to deal with the noise in the data in the absence of an underlying model.

First, instead of using the probability satisfaction as evaluator for the best formula, we design a discrimination function based on the quantitative semantics (or robustness) of STL and in particular the average robustness introduced in~\cite{BartocciBNS15}.
The average robustness enables us to differentiate among STL classifiers that have similar discriminating performance with respect to the data by choosing the most robust one.  This gives us more information than just having the probability of satisfaction: for each trajectory, we can evaluate how much is it closed to the  ``boundary''  of the classifier, instead of only checking whether it satisfies or not a formula.
We then modify the discrimination function and the GP-UCB algorithm used in~\cite{BBS14} and~\cite{BufoBSBLB14} to better deal with noisy data and to use the quantitative semantics to emulate the average robustness distribution with respect to the parameter space of the formula.

Second, we reduce the computational cost of the \emph{Evolutionary Algorithm} (EA) by using a lightweight configuration (i.e., a low threshold of max number of iterations) of the GP-UCB optimization algorithm to estimate the parameters of the formulas at each generation.
The EA algorithm generates, as a final result, an STL formula tailored for classification purpose.

We compare our framework with our previous methodology
\cite{BufoBSBLB14} and with the decision-tree based approach presented in~\cite{bombara_decision_2016} on an anomalous trajectory detection problem of a naval surveillance and on an Assisted Ventilation in Intensive Care Patients.
Our experiments indicate that the proposed approach
outperforms our previous work with respect to accuracy and show that it produces in general more compact, and easy to read, temporal logic specifications with respect to the decision-tree based  approach with a comparable speed and accuracy.

\paragraph*{\textbf{Paper structure}}
The rest of the paper is organized as follows: in the next section we present the {\it Signal Temporal Logic} and its robust semantics. We then introduce the problem considered in Section  \ref{sec:problem}. In section \ref{sec:methodology}, we describe our approach. The results are presented in Section \ref{sec:results}. Finally, we conclude the paper in Section \ref{sec:conclusion}, by discussing the implications of our contribution, both from the practical and the methodological aspect and some directions of improvement.
\section{Methodology}
\label{sec:methodology}

Learning an STL formula can be separated in two optimization problems: the learning of the formula structure and the synthesis of the formula parameters.
The structural learning is treated as a discrete optimization problem using an {\it Genetic Algorithm} (GA); the parameter synthesis, instead, considers a continuous parameter space and exploits an active learning algorithm, called {\it Gaussian Process Upper Confidence Bound (GP-UCB)}.
The two techniques are not used separately but combined together in a bi-level optimization. The GA acts externally by defining a set of template formulae. Then, the GP-UCB algorithm, which acts at the inner level, finds the best parameter configuration such that each template formula better discriminates between the two datasets.
For both optimizations, we need to define a score function to optimize, encoding the criterion to discriminate between the two datasets.

\begin{algorithm}[htbp]
\caption{ROGE -- RObustness GEnetic }
\label {algo:ROGE}
\vspace{1mm}
\begin{algorithmic}[1]
\REQUIRE $\mathcal{D}_p, \mathcal{D}_n, \mathbb{K}, Ne, Ng, \alpha, s$
\STATE $gen \gets \textsc{generateInitialFormulae}(Ne,s)$
\STATE $gen_{\Theta} \gets \textsc{learningParameters}(gen, G, \mathbb{K})$
\FOR {$i=1 \dots Ng$}
	\STATE $subg_{\Theta} \gets \textsc{sample}(gen_{\Theta} ,F)$
	\STATE $newg \gets \textsc{evolve}(subg_{\Theta},\alpha)$
    \STATE $newg_{\Theta} \gets \textsc{learningParameters}(newg, G, \mathbb{K})$
	\STATE $gen_{\Theta} \gets \textsc{sample}(newg_{\Theta} \cup gen_{\Theta}, F)$
\ENDFOR
\RETURN $gen_{\Theta}$
\end{algorithmic}
\end{algorithm}

Our implementation, called \emph{RObustness GEnetic} (ROGE) algorithm is described in  Algorithm \ref{algo:ROGE}.  First, we give an overview of it and then we described each specific function in the next subsections.
The algorithm requires as input the dataset $\mathcal{D}_p$(good) and $\mathcal{D}_n$(bad), the parameter space $\mathbb{K}$, with the bound of each considered variable, the size of the initial set of formulae $Ne$, the number of maximum iterations $Ng$, the mutation probability $\alpha$ and the maximum initial formula size $s$.
The algorithm starts generating (line 1) an initial set of PSTL formulae, called $gen$.
For each of these formulae (line 2), the algorithm learns the best parameters accordingly to a discrimination function $G$.  We call $gen_{\Theta}$ the generation for which we know the best parameters of each formula.
During each iteration, the algorithm (line 4), guided by a fitness function F , extracts  a subset $subg_{\Theta}$ composed by the best $Ne/2$ formulae, from the initial set $gen_{\Theta}$.
From this subset (line 5), a new set $newg$ composed of $Ne$ formulae is created by employing the \textsc{Evolve} routine, which implements two {\it genetic} operators.
Then (line 6), as in line 2, the algorithm identifies the best parameters of each formula belonging to $newg$.
The new generation $newg_{\Theta}$ and the old generation $gen_{\Theta}$ are then merged together (line 7). From this set of $2Ne$ formulae the algorithm extracts, with respect to the fitness function $F$, the new generation $gen_{\Theta}$ of $Ne$ formulae.
At the end of the iterations or when the  $\textsc{stop}$ criterion is true (lines 8-12), the algorithm returns the last generated formulae.
The best formula is the one with the highest value of the discrimination function, i.e., $\phi_{best} =  \mathrm{arg max}_{\phi_{ \boldsymbol{\theta}} \in gen_{\Theta}}(G(\phi_{\boldsymbol{\theta}}))$.
We describe below in order: the discrimination function algorithm, the learning of the parameters of a formula template and the learning of the formula structure.

\subsection{\bf Discrimination Function $G(\phi)$}
We have two datasets $\mathcal{D}_p$ and $\mathcal{D}_n$ and we search for the formula $\phi$  that better separates these two classes. We define a function able to discriminate between this two datasets, such that maximising this \emph{discrimination function} corresponds to finding the best formula classifier.
In this paper, we decide to use, as evaluation of satisfaction of each formula, the quantitative semantics of STL.
As described in Section~\ref{sec:background}, this semantics computes a real-value (robustness degree) instead of only a yes/no answer.

Given a  dataset $\mathcal{D}_i$, we  assume that the data comes from an unknown stochastic process $\vec X_i$. The process in this case is like a black-box for which we observe only a subset of trajectories, the dataset $\mathcal{D}_i$.
We can then estimate the averages robustness $\bb{E}(R_{\phi}|\vec X_i) $ and its variance $\sigma^2(R_{\phi} | \vec X_i)$, averaging over $\mathcal{D}_i$.

To discriminate between the two dataset $\mathcal{D}_p$ and $\mathcal{D}_n$, we search the formula $\phi$ that maximizes the difference between the average robustness of  $\vec X_p$, $\bb{E}(R_{\phi} | \vec X_p) $, and the average robustness of  $\vec X_n$, $\bb{E}(R_{\phi}| \vec X_n) $ divided by the sum of the respective standard deviation
\begin{equation}
\label{eqn:discrfcn}
G({\phi}) =  \frac{\bb{E}(R_{\phi} | \vec X_p) - \bb{E}(R_{\phi} | \vec X_n)}{\sigma(R_{\phi} | \vec X_p) + \sigma(R_{\phi} | \vec X_n) }.
\end{equation}
This formula is proportional to the probability that a new point sampled from the distribution generating $\vec X_p$ or  $\vec X_n$, will belong to one set and not to the other. In fact, an higher value of $G({\phi})$ implies that the two average robustness will be sufficiently distant relative to their  length-scale, given by the standard deviation $\sigma$.

As said above, we can  evaluate only a statistical approximation of $G({\phi})$ because we have a limited number of samples belonging to $\vec X_p$ and  $\vec X_n$.
We will see in the next paragraph how we tackle this problem.

\subsection{\bf GB-UCP: learning the parameters of the formula}
\label{sec:methodology:GPUCB}
Given a formula $\phi$ and a parameter space  $ \mathbb{K}$, we want to find the parameter configuration $\boldsymbol{\theta} \in  \mathbb{K}$ that maximises the score function $g (\boldsymbol{\theta}) := G(\phi_{\boldsymbol{\theta}})$, considering that the evaluations of this function are noisy. So, the question is how to best estimate (and optimize) an unknown function from observations of its value at a finite set of input points. This is a classic non-linear non-convex optimization problem that we tackle by means of
the GP-UCB algorithm \cite{gpucb}.
This algorithm interpolates the noisy observations using a stochastic process
(a procedure called emulation in statistics) and optimize the emulated function using the uncertainty of the fit to determine regions where the true maximum can lie.
More specifically, the GP-UCB  bases its emulation phase  on Gaussian Processes, a Bayesian non-parametric regression approach~\cite{Rasmussen2007}, adding candidate maximum points to the training set of the GP in an iterative fashion, and terminating when no improvement is possible (see \cite{gpucb} for more details).

After this optimization, we have found a formula that separates the two datasets, not from the point of view of the satisfaction (yes/no) of the formula but from the point of view of the robustness value. In other words, there is a threshold value $\alpha$ such that  $\bb{E}(R_{\phi}|\vec X_p) > \alpha$ and $\bb{E}(R_{\phi}|\vec X_n) \leq \alpha$. Now, we consider the new STL formula $\phi'$ obtained translating the atomic predicates of $\phi$ by $\alpha$ (e.g., $y(x)>0$ becomes $y(x)-\alpha>0$). Taking into account that the quantitative robustness is achieved by combination of $\min$, $\max$, $\inf$ and $\sup$, which are linear algebraic operators with respect to translations (e.g, $\min(x,y)\pm c = \min(x \pm c,y \pm c)$), we easily obtain that $\bb{E}(R_{\phi'}|\vec X_p) > 0$ and $\bb{E}(R_{\phi'}|\vec X_b)< 0$. Therefore, $\phi'$ will divide the two datasets also from the point of view of the satisfaction.

\subsection{\bf Genetic Algorithm: learning the structure of the formula}
 To learn the formula structure, we exploit a modified version of the Genetic Algorithm (GA) presented in  \cite{BufoBSBLB14}. GAs belongs to the larger class of evolutionary algorithms (EA). They are used for search and optimization problems. The strategy takes inspiration from the genetic area, in particular in the selection strategy of species.  Let us see now in detail the genetic routines of the ROGE algorithm.

\paragraph*{\bf $gen = \textsc{generateInitialFormulae}(Ne,s)$.} This routine generates the initial set of $Ne$ formulae. A fraction $n_l < Ne$ of this initial set is constructed by a subset of the temporal properties: $\ev{I}\mu$,  $\glob{I}\mu$, $\mu_1\until{I}\mu_2$, where the atomic predicates are of the form $\mu = (x_i > k)$  or $\mu = (x_i \le k)$ or simple boolean combinations of them (e.g. $\mu = (x_i > k_i) \wedge (x_j > k_j)$). The size of this initial set is exponential accordingly to the number of considered variables $x_i$. If we have few variables we can keep all the generated formulae, whereas if the number of variables is large we consider only a random subset. The remain $n_r = Ne - n_l$ formulae are chosen randomly from the set of formulae with a maximum size defined by the input parameter $s$. This size can be adapted to have a wider range of initial formulae.

\paragraph*{\bf $ subg_{\Theta} = \textsc{sample}(gen_{\Theta} ,F)$.}
This procedure extracts from $gen_{\Theta}$ a subset $subg_{\Theta}$ of $Ne/2$ formulae, according to a fitness function $F(\phi) = G(\phi) - S(\phi)$.
The first factor $G(\phi)$ is the dicrimination function previously defined and $S(\phi)$ is a size penalty, i.e. a function penalizes formulae with large sizes.  In this paper, we consider $S(\phi) = g \cdot p^{size(\phi)}$, where $p$ is heuristically set such that $p^5 = 0.5$, i.e. formulae of size 5 get a 50\% penalty, and $g$ is adaptively computed as the average discrimination in the current generation. An alternative choice of $p$ can be done by cross-validation.

\paragraph*{\bf $ newg = \textsc{evolve}(subg_{\Theta},\alpha)$.}
This routine defines a new set of formulae ($newg$)  starting from $subg_{\Theta}$, exploiting  two {\it genetic} operators: the {\it recombination}  and the {\it mutation} operator.
The recombination operator takes as input two formulae $\phi_1$, $\phi_2$ (the parents), it randomly chooses a subtree from each formula and it swaps them, i.e. it assigns the subtree of $\phi_1$ to $\phi_2$ and viceversa. As a result, the two generated formulae (the children) share different subtrees of the parents' formulae.
The mutation operator takes as input one formula and induce a random change (such as inequality flip, temporal operator substitution, etc.) on a randomly selected node of its tree-structure. The purpose of the genetic operators is to introduce innovation in the offspring population which leads to the optimization of a target function ($G$ in this case). In particular,  recombination is related to exploitation, whereas  mutation to exploration.
The \textsc{evolve} routine implements an iterative loop that at each iteration  selects which genetic operators to apply. A number $p \in [0,1]$ is randomly sampled.  If its value is lower than the mutation probability, i.e.,  $p \le \alpha$, then the mutation is selected, otherwise a recombination is performed.
At this point the input formulae of the selected genetic operator are chosen randomly between the formulas composing $subg_{\Theta}$ and the genetic operators are applied.
In our implementation, we give more importance to the exploitation; therefore the mutation acts less frequently than the recombination (i.e., $\alpha \le 0.1$).
The iteration loops will stop when the number of generated formula is equal to $Ne$, i.e. twice the cardinality of $subg_\Theta$.

\section{Problem formulation}
\label{sec:problem}

In this paper, we focus our attention on learning the best property
(or set of properties) that  discriminates trajectories belonging to
two different classes, say good and bad, starting from a labeled
dataset of observed trajectories. Essentially, we want to tackle
a supervised two-class classification problem over trajectories,
by learning a temporal logic discriminant, describing the temporal
patterns that better separate two sets of observed trajectories.

The idea behind this approach is that there exists an unknown
procedure that, given a temporal trajectory, is able to decide if
the signal is good or bad. This procedure can correspond to
many different things, e.g., to the reason of the failure of
a sensor that breaks when it receives certain inputs.
In general, as there may not be an STL specification
 that perfectly explains/mimics the unknown procedure,
our task is to approximate it with the most effective one.

Our approach works directly with observed data, and
avoids the reconstruction of an intermediate generative
model $p(\vec x|z)$ of trajectories $x$ conditioned on
their class $z$, as in \cite{BufoBSBLB14,BBS14}.
The reason is that such models can be hard to construct,
even if they provide a powerful regularization, as they
enable the generation of an arbitrary number of samples
to train the logic classifier.

In a purely data-driven setting, to build an effective classifier,
we need to consider that training data is available in limited
amounts and it is typically noisy.  This reflects in the necessity
of finding methods that guarantee good generalization
performance and avoid overfitting.
In our context, overfitting can result in overly complex formulae,
de facto encoding the training set itself rather than guessing
the underlying patterns that separate the trajectories.  This can
be faced by using a score function based on robustness of
temporal properties, combined with suitably designed
regularizing terms.

We want to stress that the approach we present here,
due to the use of the average robustness of STL properties,
can be easily tailored to different problems,  like finding the
property that best characterise a single set of observations.

\section{RELATED WORK}

Mining temporal logic requirements is an emerging field of research in the analysis
of cyber-physical systems (CPS)~\cite{Xu2018,AckermannCHRSL10,AsarinDMN11,HoxhaDF18,BBS14,JinDDS15,nguyen_abnormal_2017,ZhouRWT17,kong_temporal_2017,bombara_decision_2016,BufoBSBLB14}.
This approach is orthogonal to
active automata learning (AAL) such as $L*$ Angluin's
algorithm~\cite{Angluin87} and its recent variants~\cite{IsbernerHS14,SteffenHI12}.
AAL is suitable to capture the behaviours of black-box reactive systems
and it has been successfully  employed in the field of
CPS to learn how to interact with the
surrounding environments~\cite{Chen2013,FuTHC14}.  Mining temporal logic
requirements has the following advantages with respect to AAL.
The first is that it does not require to interact with a reactive system.
AAL needs to query the system in order to learn a Mealy machine
representing the relation between the input provided and the output observed.
Mining temporal logic requirements can
be applied directly to a set of observed signals without
the necessity to provide an input.
The second is the possibility to use temporal logic
requirements within popular tools (such as Breach~\cite{Donze10}
and S-TaLiRo~\cite{AnnpureddyLFS11}) for monitoring and falsification
analysis of CPS models.

Most of the literature related to  temporal logic inference from data
focuses in particular on the problem of learning the optimal
parameters given a specific template formula~\cite{Xu2018,AsarinDMN11,HoxhaDF18,BBS14,JinDDS15,nguyen_abnormal_2017,ZhouRWT17}.
In~\cite{AsarinDMN11}, Asarin et al. extend {\it the Signal Temporal Logic} (STL)~\cite{Maler2004}
with the possibility to express the time bounds of
the temporal operators and the constants of the
inequalities as parameters. They also
provide a geometric approach to identify the subset of the
parameter space that makes a particular signal
to satisfy an STL specification.
 Xu et al. have recently proposed in~\cite{Xu2018}  a temporal
logic framework called \emph{CensusSTL}
for multi-agent systems that consists of
\emph{an outer logic STL formula}  with a variable
in the predicate representing the number of
agents satisfying \emph{an inner
logic STL formula}.  In the same paper the
authors propose also a new inference algorithm
similar to~\cite{AsarinDMN11} that given the
templates for both the \emph{inner} and \emph{outer}
formulas,  searches for the optimal parameter values
that make the two formulas capturing the trajectory data
of a group of agents.
In~\cite{HoxhaDF18} the authors use the same parametric STL extension in
combination with the quantitative semantics of STL to
perform a counter-example guided inductive parameter
synthesis.  This approach consists in iteratively
generating a counterexample by executing a
falsification tool for a template formula.
The counterexample found at each step is then
used to further refine the parameter set and
the procedure terminates when no other
counterexamples are found.
In general, all these methods, when working directly
with raw data, are potentially vulnerable to the
noise of the measurements and they are limited
by the amount of data available.

Learning both the structure and the parameters of a
formula from a dataset poses even more challenges~\cite{kong_temporal_2017,bombara_decision_2016,BBS14,BufoBSBLB14}.
This problem is usually addressed in two steps,
learning the structure of the formula and synthesizing its parameters.
In particular, in~\cite{kong_temporal_2017} the structure
of the formula is learned by exploring a directed acyclic graph and the method  exploits {\it Support Vector Machine} (SVM)
for the parameter optimization.  In~\cite{bombara_decision_2016} the authors use instead a \emph{decision tree} based approach for learning the formula, while the optimality is evaluated using heuristic
impurity measures.

In our previous works~\cite{BBS14,BufoBSBLB14} we have also
addressed the problem of learning both the structure and the
parameters  of a temporal logic specification from data.  In~\cite{BBS14} the structure
of the formula is learned using a heuristics algorithm, while
in~\cite{BufoBSBLB14} using a genetic algorithm. The
synthesis of the parameters is instead performed in both cases
exploiting the {\it Gaussian Process Upper Confidence Bound}
(GP-UCB)~\cite{gpucb} algorithm, statistically emulating
the satisfaction probability of a formula for a given set of parameters.
In both these methodologies, it is required to learn first a statistical
model from the training set of trajectories. However, the statistical
learning of this model can be a very difficult task and this is one of
the main reason for proposing our new approach.

\section{Case studies and Experimental Results}
\label{sec:results}

\subsection{Maritime Surveillance}
As first case study, We consider the maritime surveillance problem presented in~\cite{bombara_decision_2016} to compare our framework with their {\it Decision Tree} (DTL4STL) approach.  The experiments with the DTL4STL approach were implemented
in \texttt{Matlab}, the code is available at~\cite{DTL4STL}.
We also test our previous implementation presented
in~\cite{BufoBSBLB14} with the same benchmark.
Both the new an the previous learning procedure
were implemented in \texttt{Java} (JDK 1.8\_0) and run on a Dell XPS, Windows 10 Pro, Intel Core i7-7700HQ 2.8 GHz,
8GB 1600 MHz memory.

The synthetic dataset of naval surveillance reported
in~\cite{bombara_decision_2016} consists of 2-dimensional
coordinates traces of vessels behaviours.  It considers two kind
of anomalous trajectories and regular trajectories, as illustrated in Fig.~\ref{fig:naval dataset}.  The dataset contains 2000 total trajectories (1000 normal and 1000 anomalous) with 61 sample points per trace.
   \begin{figure}[thpb]
      \centering

       \includegraphics[width=1.1 \textwidth]{navalPlot.eps}

      \caption{The regular trajectories are represented in green. The anomalous trajectories which are of two kinds, are represented respectively in blue and red.}
      \label{fig:naval dataset}
   \end{figure}
We run the experiments using a 10-fold cross-validation in order to collect the mean and variance of the misclassified trajectories of the validation set.
Results are summarized in Table~\ref{tab:navalResults}, where we report also the average execution time.
We test also our previous implementation~\cite{BufoBSBLB14}
without a generative model from data. It performs so
poorly on the chosen benchmark that is not meaningful to report it: the misclassification rate on the validation set is around 50\%.
In Table \ref{tab:navalResults}, we also report DTL4STL$_p$ the DTL4STL performance declared in~\cite{bombara_decision_2016}, but we were not be able to reproduce them in our setting.

\begin{table}[!t]
\begin{center}
\scalebox{1.15}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
& ROGE &DTL4STL & DTL4STL$_p$ \\
\hline
Mis. Rate &$0$  &  $0.01 \pm 0.013$& $0.007\pm0.008$ \\
\hline
Comp. Time &$73 \pm18$ & $144 \pm24$& -  \\

\hline
\end{tabular}}
\end{center}
\caption{The comparison of the computational time (in sec.), the mean missclassification rate and its standard deviation between the learning procedure using the RObust GEnetic algorithm, the Decision-Tree (DTL4STL) Matlab code provided by the authors and the results reported in  \cite{bombara_decision_2016} (DTL4STL$_p$).}
\label{tab:navalResults}
\end{table}
In terms of accuracy our approach is comparable with respect to the performance of the DTL4STL. In terms of computational cost, the decision tree approach is slightly slower but it is implemented in \texttt{Matlab} rather than \texttt{Java}.

An example of formula that we learn with ROGE is
\begin{equation}
\label{eqn:bestourformula}
\phi=((x_2 > 22.46)\, \mathcal{U}_{[49, 287]}\, (x_1 \le 31.65))
\end{equation}
The formula (\ref{eqn:bestourformula}) identifies all the regular trajectories, remarkably resulting in a misclassification error equal to zero, as reported in Table \ref{tab:navalResults}. The red anomalous trajectories falsify the predicate $x_2 > 22.46$ before verifying $x_1\le31.65$, on the contrary the blue anomalous trajectories globally satisfy $x_2 > 22.46$ but never verify $x_1 \le 31.65$ (consider that all the vessels start from the top right part of the graphic in Figure \ref{fig:naval dataset}). Both these conditions result in the falsification of Formula(~\ref{eqn:bestourformula}), which on the contrary is satisfied by all the regular trajectories. In Figure \ref{fig:naval dataset}, we have reported the threshold $x_2=22.46$ and $x_1=31.65$.

An example instead of formula found by the DTL4STL tool using the same dataset is the following:
\begin{eqnarray*}
\small
\label{eqn:bestDTformula}
\psi =
(((\mathcal{G}_{[187,196)}x_{1}<19.8) \wedge
(\mathcal{F}_{[55.3,298)}x_{1}>40.8) ) \vee ((\mathcal{F}_{[187,196)}x_{1}>19.8) \wedge \\ ((\mathcal{G}_{[94.9,296)}x_{2}<32.2)  \vee ((\mathcal{F}_{[94.9,296)}x_{2}>32.2) \wedge (((\mathcal{G}_{[50.2,274)}x_{2}>29.6) \wedge\\
  (\mathcal{G}_{[125,222)}x_{1}<47) ) \vee
   ((\mathcal{F}_{[50.2,274)}x_{2}<29.6) \wedge (\mathcal{G}_{[206,233)}x_{1}<16.7) )))))
\end{eqnarray*}

The specific formula generated using ROGE is simpler than the formula generated using DTL4STL and indeed it is easier to understand it. Furthermore DTL4STL does not consider the until operator in the set of primitives (see \cite{bombara_decision_2016} for more details).

\subsection{Ineffective Inspiratory Effort (IIE)}
The Ineffective Inspiratory Effort (IIE) is one of the major problems concerning the mechanical ventilation of patients in intensive care suffering from respiratory failure.
The mechanical ventilation uses a pump to help the patient in the process of ispiration and expiration, controlling the $\mathit{flow}$
of air into the lugs and the airway pressure ($paw$).
Inspiration is characterised by growth of pressure up to the selected $paw$ value and positive $\mathit{flow}$, while expiration has a negative $\mathit{flow}$ and a drop $paw$. The exact dynamics of these  respiratory curves stricly depends on patient and on ventilatory modes. We can see an example in~Fig.~\ref{fig:breaths} of two regular (blue regions) and one ineffective (red region) breaths. An IIE occurs when the patients tries to inspire, but its effort is not sufficient to trigger the appropriate reaction of the ventilator. This results in a longer expiration phase. Persistence of IIE may cause permanent damages of the respiratory system.

An early detection of anomalies using low-cost methods is a key challenge to prevent IIE conditions and still an open problem.
\begin{figure}[thpb]
      \centering
      \includegraphics[width=1.\textwidth]{breaths.png}
      \caption{Example of regular (blue regions) and ineffective (red region) breaths.}
      \label{fig:breaths}
  \end{figure}
In~\cite{BufoBSBLB14}, starting with a dataset of discrete time series and sampled $\mathit{flow}$ values (labeled good and bad) from a single patient, the authors statically designed generative models of $\mathit{flow}$ and of its derivative $\mathit{flow}'$, for regular and ineffective signals. Then they used the simulations of such models to identify the best formula/formulae discriminating between them. In that contribution trajectories with different lengths are considered, treating as false trajectories that are too short to verify a given formula, and use this to detect the length of trajectories and separate anomalous breaths which last longer than regular ones.
However, using the information about the duration of a breath is not advisable if the goal is to monitor the signals at runtime and detect the IEE at their onset.
For this reason, in our approach we consider a new dataset, cutting the original trajectories used to generate the stochastic model in Bufo et al.~\cite{BufoBSBLB14} to a maximum common length of the order of 2 seconds. We also apply a moving average filter to reduce the noise in the computation of $\mathit{flow}'$.

Specifically, the new dataset consists of 2-dimensional traces of $\mathit{flow}$ and its derivative, $\mathit{flow}'$, containing a total of 288 trajectories (144 normal and 144 anomalous) with 134 sample points per trace. The time scale is in hundredths of a second.
We run the experiments using a 6-fold cross-validation and report our results and comparison on the new dataset with DTL4STL~\cite{bombara_decision_2016} in Table~\ref{tab:breathResults}.

\begin{table}[!t]
\begin{center}
\scalebox{1.15}{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
& ROGE & DTL4STL \\
\hline
Mis. Rate & $0.17\pm 0.01$  & $0.23\pm0.07$  \\
\hline
False. Pos & $0.20\pm 0.02$  & $0.23\pm0.07$ \\
\hline
False. Neg & $0.14\pm 0.02$  & $0.20\pm0.15$ \\
\hline
Comp. Time & $65\pm7$ & $201 \pm 7$  \\
\hline
\end{tabular}}
\end{center}
\caption{The comparison of the computational time (in sec.), the mean missclassification rate and its standard deviation between the learning procedure using the RObust GEnetic algorithm, }
\label{tab:breathResults}
\end{table}

An example of formula that we learn with ROGE is:
\begin{equation}
\label{eqn:bestourformulaIIE1}
\phi =  \ev{[31, 130]}( (\mathit{flow} \ge -670) \vee (\mathit{flow}' \le -94))
\end{equation}

Formula $\phi$ identifies the anomalous trajectories, stating that at a time between $31$ sec and $130$ seconds either $\mathit{flow}$ is higher than $-670$ or  $\mathit{flow}'$ is below $-94$.
This is in accordance with the informal description of an IIE, principally caused by an unexpected increment of the $\mathit{flow}$ during the expiration phase followed by a rapid decrease. Indeed, one of the main characteristic of an IIE is the presence of a small hump in the $\mathit{flow}$ curve during this phase. In Fig.~\ref{fig:hump} we report some of the trajectories of the dataset, showing the areas where the property is satisfied.
     \begin{figure}[thpb]
      \centering      \includegraphics[width=1.\textwidth]{breathprop.png}

      \caption{  100 regular (green) and  100 ineffective (red) flow  and $\mathit{flow}'$ trajectories during the expiration phase. The light blue rectangles correspond to the satisfaction area of formula~(\ref{eqn:bestourformulaIIE1}). One regular (green) and one ineffective (red) trajectories are showed thicker.}
      \label{fig:hump}
   \end{figure}

On average, formulae found by the DTL4STL tool on the new dataset are disjunctions or conjunctions of$~10$ eventually or always subformulae, which are not readily interpretable.

Our approach on the new dataset is better in term of accuracy and computation time with an improvement of 22\% and 67\%, respectively.

Similarly to the previous test case, we compare our approach with \cite{BufoBSBLB14}, performed directly on the dataset, and not on the generative model. That approach performs so poorly also in this benchmark, obtaining a misclassification rate of $0.47$, which is comparable with a random classifier. The problem here is that the methods proposed in \cite{BufoBSBLB14}, differently from the one presented here, relies on a large number of model simulations, and it is not suited to work directly with observed data.

If we give up to online monitoring and consider only full breaths, we can improve classification by rescaling their duration into [0,1], so that each breath lasts the same amount of time. In this case, we learn a formula with misclassification rate of $0.05\pm0.01$, while DTL4STL reaches a misclassification rate of $0.07\pm0.02$. This suggests that the high variability of durations of ineffective breaths has to be treated more carefully.

\section*{Acknowledgment}

E.B.\ and L.N.\ acknowledge the partial support of the Austrian National
Research Network  S 11405-N23 (RiSE/SHiNE) of the Austrian Science
Fund (FWF). E.B.,\ L.N.\ and S.S.\ acknowledge the partial support of the ICT COST Action IC1402 (ARVI).

\bibliographystyle{splncs}