\section{Introduction}
\label{sec:Introduction}

\subsection{Aperture synthesis}
The Atacama Large Millimeter/submillimeter Array (ALMA) is revolutionizing  observational astrophysics. With its 66 antennas located on the Atacama desert it provided the sharpest ever images of the submillimeter sky, for example, images of the protoplanetary disks at 1 au resolution \cite{Andrews.Wilner.ea2016}. In order to obtain such a resolution at a distance to a nearby star-forming region at 1.3 mm a telescope diameter of $\sim 15$ km are needed. Since the construction challenges of such an antenna, especially if one would like to make it steerable, are far beyond current technical capabilities, in radio astronomy domain we often turn to aperture synthesis techniques, where instead of a single dish, a combination of smaller antennas is used, and with interference of signal between each antennas a resolution compared to the a telescope of a size of the greatest distance between the two antennas in an array (i.e. baseline) is achieved.

This does not come without a cost: sampling the baselines is never complete compared with a single dish telescopes. This means we do not receive complete information at all baselines and therefore to create an image of the sky we are operating with missing information.

A direct measurement of an interferometer is the interference pattern between two given antennas. This pattern is related to the sky brightness observed by the antennas. The recorded complex value called {\it visibility} is then a Fourier transform of the sky brightness, the quantity which observations aim to recover. Therefore, a simplified imaging process consists of a (reverse) Fourier transformation of the measured visibilities (while the non-measured visibilities are set to 0) to obtain first approximation of the sky brightness, the so-called {\it dirty image} \citep{Jackson2008,Thompson.Moran.ea2001}. Once this is achieved, it becomes apparent that even with modern interferometer like ALMA with many baselines sampling the so-called UV plane we achieve a rather poor quality image. Refining this image is the topic of this work.

\subsection{CLEAN as fa standard approach to the imaging of the interferometric data}

A standard technique to improve the image quality of the interferometric observations has been developed by Hogbom in 1974 \cite{Hoegbom1974} and is called \texttt{clean}. \texttt{Clean} makes use of the well-defined point-spread function of a given antenna configuration. The algorithm identifies the point sources in the initial dirty image, point sources are then approximated with a delta function, convolved it with a dirty beam (i.e. the assumed pattern that the point-source would create on the dirty image), scales it with brightness of the suspected point source and subtracts the point-source patter from the dirty image. This is an iterative process which ideally ends when all point sources are removed from the dirty image so that it consists only of noise. \citep{Jackson2008}

One issue of \texttt{clean} that can be easily identified is in the assumption that the sky is composed of point sources. This results in \texttt{clean} struggling with imaging of extended sky brightness structures. Several modifications to the original \texttt{clean} algorithm have been implemented to mitigate this issue, such as multi-scale clean which allows to set a point-source that would be a Gaussian rather than a delta function \cite{Junklewitz.Bell.ea2015,Cornwell2008}.

There are several steps in the process of cleaning which can be modified to improve the process. Masking is a method to restrict the area where the algorithm will look for point sources to a selected region on the sky. Weighting allows to attribute different weights to different u,v scales, allowing for a trade-off between resolution and sensitivity.

Although \texttt{clean} became a gold standard of the interferometric data imaging, producing many stunning images of the astronomical objects its limitation beg to seek alternatives, especially in cases where its assumptions are not met.

\subsection{\texttt{resolve} algorithm and IFT}
Imaging of the interferometric data can be presented as an inference problem, since we operate on an incomplete measurement problem, trying to find the true sky emission from the received data. Radio Extended SOurces Lognormal deconVolution Estimator \texttt{resolve} \footnote{https://gitlab.mpcdf.mpg.de/ift/resolve} \cite{Junklewitz.Bell.ea2015,Arras.Knollmueller.ea2018} is designed in the Information Field Theory (IFT) framework \citep{Ensslin.Frommert.ea2009}. IFT enables to use Bayesian inference methods in the context of mathematical framework of field theory. This is well fitted to the issue of imaging the sky brightness. IFT algorithms are implemented in \texttt{resolve} through Python package NIFTy \footnote{https://gitlab.mpcdf.mpg.de/ift/NIFTy} \cite{Selig.Bell.ea2013,Arras.Baltac.ea2019}.

The measurement equation can be presented as: $d = R(s) + n$, where $R$ is a response of the instrument to the original physical signal ($s$), and $n$ is noise. In the IFT framework, obtained data $d$ is analyzed in order to find the most probable form of $s$, which is signal from the field, in this case, the sky. This takes the form of Bayes' theorem:
\begin{equation}
     P(s|d) =\frac{P(d|s)P(s)}{P(d)},
\end{equation}
$P(d|s)$ - likelihood that a given data (d) has been produced from a signal (s), $P(s)$ - is a prior knowledge about the signal, and P(d) is a normalization factor.

\texttt{Resolve} is developed in order to optimize imaging of extended and diffuse radio sources and to provide reliable noise estimation \cite{Junklewitz.Bell.ea2015}. It has been successfully used for example in imaging of Cygnus A \citep{Arras.Bester.ea2021} observed with the VLA and the M87* black hole environment with VLBI \citep{Arras.Frank.ea2022}. In this work we want to explore \texttt{resolve} capabilities for ALMA interferometer, specifically compared with the most commonly used \texttt{clean} algorithm.

\section{Test datasets}
In this section we present data used to test interferometric imaging methods: \texttt{resolve} and \texttt{clean}. We use a simple simulated dataset in order to have a  complete control over the input information, as well as the real ALMA observations.
\subsection{Simulated data}
\label{sec:simulated_data}

In this section, we describe the preparation of the simulated data. The major advantage in using simulated datasets to analyze the imaging procedures is a complete control of the input parameters. We control perfectly the shape and brightness of the sources, as well as configuration and behaviour of the telescope array. On the downside of this approach is the difficulty to realistically model the noise acquired during observation and simplistic assumptions about the sky brightness.

We generate a simple 2D array consisting of five Gaussian components of different brightness, size and position angle. This model is input into the Common Astronomy Software Applications package (CASA; \cite{McMullin.Waters.ea2007}) task \texttt{simalma}. This task converts image to the sky model, imposing physical properties to the sky, such as spherical coordinates, pixel dimensions, field-of-view and brightness in physical units. Fig. \ref{fig:simalma_skymodel} (left) shows the sky model generated with \texttt{simalma} task in CASA.
Afterwards, the task is simulating observations of the given sky model with ALMA observatory. For this specific case we simulate observations with ALMA configuration C-3 at 230 GHz (ALMA band 6), which results in effective resolution of 0.7", determined approximately by largest available baseline (distance between two antennas). In the case of C-3 configuration this is $\sim$ 500 m.

The \texttt{simalma} task returns a calibrated Measurement Set (MS), that consists of complex visibilities. Those visibilities are Fourier transform of the sky brightness, therefore reverse Fourier transformation provides a dirty image of the observed sky. Further on we describe the process of imaging with \texttt{resolve} and \texttt{tclean}.

\begin{figure}

    \includegraphics[width=0.5\textwidth]{sim314_skymodel.png}
    \includegraphics[width=0.45\textwidth]{Sz114_best_v3.png}
    \caption{{\it Left:} Sky model which serves as an input for creating simulated ALMA Measurement Set. {\it Right:} Best \texttt{tclean} image of Sz114 protoplanetary disk from ALMA observations.}
    \label{fig:simalma_skymodel}
\end{figure}

\subsection{ALMA data}

Here we describe archival ALMA observations used to test application of \texttt{resolve} on real-life example. It is important to compare it with a well-understood case, in which the calibration and imaging examples with other tools are already available. We select a protoplanetary disk Sz114 observed within DSHARP ALMA Large Program \citep{Andrews.Huang.ea2018a}. This was a milestone program for ALMA observatory, showcasing the richness of substructures within disks, often associated with ongoing planet formation. It is therefore especially interesting to search for innovative methods of imaging those data in order to verify current conclusions, as well as open avenue for new studies.

The Sz114 disk shows relatively smooth structure compared with other extreme cases and therefore it serves well as a test case to try to identify underlying structure with \texttt{resolve}. Fig. \ref{fig:simalma_skymodel} (right) presents image of the disk obtained with \texttt{tclean} as delivered by the Large Program team \citep{Andrews.Huang.ea2018a}.

\subsection{Imaging of the simulated dataset}

First, in order to image the simulated measurement set, we use \texttt{tclean} task in CASA software which implements CLEAN algorithm as in Hogbom 1974 \cite{Hoegbom1974}. We ran \texttt{tclean} without any constraint on where to look for point sources in the image (i.e. without any masking), for 20000 iterations, or until threshold of 0.3 mJy/beam was reached. Pixel size of the reconstructed image is set to 0.1" and image size is set to 512$\times$512 pixels.

This threshold is set based on the earlier, quick \texttt{tclean}, so that the algorithm does not attempt to find sources from the residual image consisting purely of noise. In the default settings, weighting of baselines is set to 'natural', which means it associate the baseline with weight proportional to the sampling density (i.e. the most covered baselines have the highest weight). Since there is much more baselines sampling the larger scales, this results in putting more weight to lower resolution, which results in achieving lower resolution than the sky model image. Therefore we also attempt an imaging with Briggs weighting with robust parameter 0.5, which moves the balance of weighting toward longer baselines increasing the resolution but decreasing signal-to-noise ratio. The result of the  \texttt{tclean} run is presented in Fig. \ref{fig:simalma_res_cl} (bottom).

For \texttt{resolve} image of the simulated data we run 30 iterations, using Newton optimizer. We did not make any assumptions on the presence of the point sources in the data. Parameters of the \texttt{resolve} run are summarized in Table \ref{tab:tab_resolve}. For detailed explanation of the parameters see Arras et al. \cite{Arras.Bester.ea2021}.
\vspace{1cm}

\begin{table}
\caption{Parameters of the \texttt{resolve} run on simulated dataset (left) and Sz114 ALMA data (right)}

\centering
\begin{minipage}{0.45\linewidth}
\label{tab:1}
\begin{tabular}{l c c c c c c  r r c c c cc}
\hline\hline
Parameter  & Mean\\
\hline
Offset &  26 \\
Zero mode  & 1$\pm0.1$ \\
Fluctuations &  5$\pm1$\\
Power spectrum slope & -2$\pm0.5$\\
Flexibility & 1.2$\pm0.4$\\
Asperity & 0.2$\pm$0.2\\

\hline

\end{tabular}
\end{minipage}
\begin{minipage}{0.45\linewidth}
\label{tab:2}
\begin{tabular}{l c c c c c c  r r c c c cc}
\hline\hline
Parameter  & Mean\\
\hline
Offset &  20 \\
Zero mode  & 1$\pm0.2$ \\
Fluctuations &  3$\pm1$\\
Power spectrum slope & -4$\pm2$\\
Flexibility & 4$\pm0.8$\\
Asperity & 2$\pm$0.8\\

\hline

\end{tabular}
\end{minipage}

    \label{tab:tab_resolve}
\end{table}

\subsection{Imaging of the real ALMA observations}

For the Sz114 protoplanetary disk observed within DSHAPR ALMA Large Program \cite{Andrews.Huang.ea2018a} we have created an image from a single spectral window, in order for a direct comparison with \texttt{resolve}, which currently does not have a capability to make images of multiple spectral windows. From the publicly available MS file we extracted spectral window 9 with \texttt{split} task in CASA and binned it into a single channel.

For this \texttt{tclean} run we implemented both standard and multiscale imaging in order to compare the outcomes, since the observed disk presents large variety of spatial scales. We create images with 0.005" and 1024$\times$1024 pixels and run 20000 iterations up to noise threshold of 0.05 mJy is reached. In case of multiscale clean we specified three scales: at 0, 7, and 28 pixels, which means the \texttt{tclean} algorithm iterate in order to find three types of sources in the data: either a point source, which corresponds to scale 0, and extended components with Gaussian shape of 7 and 28 pixels of FWHM. Results of the imaging are presented in Fig. \ref{fig:Sz114_res_cl} (bottom).

A \texttt{resolve} imaging followed the same procedure as in the simulated data, however, since the weights were obtained with from real observation, this results in more realistic treatment of the uncertainties. Parameters of this run are shown in Tab. \ref{tab:tab_resolve}.

\section{Comparison to \texttt{tclean} and \texttt{resolve}}
In this initial study of performance of the \texttt{resolve} algorithm in operating with ALMA data, we focus on key aspects of the imaging process: flux recovery -- how much flux observed on the sky is recovered by the imaging process and how accurate it is; image quality -- effective resolution and dynamic range.

\subsection{Simulations}

We compare the quality of \texttt{resolve} and \texttt{tclean} imaging of the simple simulated ALMA observation, as described in Section \ref{sec:simulated_data}. Fig. \ref{fig:Sz114_res_cl} presents \texttt{resolve} image and associated uncertainty map, and \texttt{tclean} image with two different weighting schemes: natural and Briggs weighting with robust parameter of 0.5.

We note that for \texttt{tclean} the imaging resulted with large areas with a negative flux, which do not occur in \texttt{resolve} - positive-only emission is one of the assumptions in the prior. On the uncertainty map produce by \texttt{resolve} we observe much higher confidence associated with the source positions, compared with positions where no significant emission was detected.

\begin{figure}

        \includegraphics[width=0.5\textwidth]{sim314_resolve_sci.png}
        \includegraphics[width=0.5\textwidth]{sim314_resolve_std_stdovermean.png}
                \includegraphics[width=0.5\textwidth]{sim314_tclean_natural.png}
        \includegraphics[width=0.5\textwidth]{sim314_tclean_briggs.png}
    \label{fig:simalma_res_cl}
        \caption{Comparison of the imaging of simulated ALMA observations with  \texttt{resolve}  and  \texttt{tclean}. {\it Top:} \texttt{resolve} image (left) with associated uncertainty map. {\it Bottom:} two \texttt{tclean} images, left: with natural weighting; right: with Briggs weighting and robust = 0.5.}
\end{figure}

We measure integrated fluxes of each of the Gaussian components in both imaging methods as well as in the input sky model. Result is presented in Table \ref{tab:1}. We can note that \texttt{resolve} and \texttt{tclean} both underestimate the flux of the brightest source (i.e. the Gaussian associated with the largest integrated flux). In case of \texttt{tclean} this can be attributed to imperfect recovery of the signal from the point-spread function, but it is unclear why \texttt{resolve} underpredicts the flux. More extensive runs are needed to investigate this issue.

\begin{figure}

        \includegraphics[width=0.5\textwidth]{Sz114_resolve_sci.png}
         \includegraphics[width=0.5\textwidth]{Sz114_resolve_std_stdovermean.png}
          \includegraphics[width=0.5\textwidth]{Sz114_clean_regular.png}
          \includegraphics[width=0.5\textwidth]{Sz114_clean_msc.png}
    \caption{Comparison of the imaging of Sz114 ALMA observations with resolve and tclean. {\it Top:} \texttt{resolve} image (left) with associated relative error map. {\it Bottom:} two \texttt{tclean} images, left: with standard cleaning algorithm, right: with multi-scale algorithm}
    \label{fig:Sz114_res_cl}
\end{figure}

On the other hand
\texttt{resolve} resolve reaches good accuracy on the other two measured components. We are able to provide much more reliable error bars on \texttt{resolve} by taking the standard deviation of the measured value from all the samples in the final iteration.

\begin{table}

\begin{minipage}{0.5\textwidth}
\caption{Peak flux (in Jy arcsec$^{-2}$) of different components on the simulated image.}
\begin{tabular}{l l c c c c c  r r c c c cc}
\hline\hline
Comp & model & \texttt{tclean} & \texttt{tclean}  & \texttt{resolve}\\
 & & natural  & Briggs &  &   \\
\hline

1 &  34.11 & 24.67 & 24.96 & 26.68$\pm$5.76 \\
2 &  23.31 & 21.68 & 22.19 & 23.59$\pm$6.72 \\
3 &  16.69 & 16.31 & 16.57 & 17.82$\pm$3.86\\
\hline
\label{tab:1}
\end{tabular}

\end{minipage}
\begin{minipage}{0.5\textwidth}

\caption{Peak (in Jy arcsec$^{-2}$) and integrated flux at given radii of the Sz114 disk.}
\begin{tabular}{l l c c c c c  r r c c c cc}
\hline\hline
Radius & best & \texttt {tclean} & \texttt{tclean}  & \texttt{resolve}\\
 & & hogbom  & multiscale &  &   \\
\hline
peak &  1.57 & 1.55 & 1.49 & 5.8$\pm$3.82 \\
0.06 &  8.63 & 8.85 & 8.95 & 9.48$\pm$1.07 \\
0.15 &  22.79 & 22.82 & 22.53 & 23.21$\pm$0.68 \\
0.35 &  47.45 & 47.07 & 47.24 & 47.24$\pm$2.79 \\
\hline
\label{tab:2}
\end{tabular}

\end{minipage}

\end{table}

\subsection{ALMA data}

In the case of real ALMA data it is more difficult to assess total recovered flux since we do not have information on true flux. We compare the peak and total flux integrating over different areas of the disk \ref{tab:2}. First of all, it can be noted that all \texttt{tclean} images have comparable fluxes, therefore the multiscale did not affect significantly the flux measurements.

In comparison between \texttt{resolve} and \texttt{tclean}, we can see that peak flux measured from \texttt{resolve} image is 3 times as high as the peak in \texttt{tclean} images. This can be purely gridding effect as the pixel scale of the \texttt{resolve} image is much smaller.

At the same time we note slightly comparable values of the integrated flux measured at from different radii of the disk.

At first sight it appears that \texttt{resolve} is able to create super-resolution image of the disk. We can use the uncertainty map to understand what confidence can be associated with those structures. We note on the uncertainty map that we typically achieve better than 20$\%$ confidence on the results but the map also reaches very low confidence in some areas of the map.

\section{Conclusions}

In this work we apply Bayesian inference and field theory in the framework of Information Field Theory (IFT) using the \texttt{resolve} algorithm, to create images from radio intereferometric observations obtained with ALMA array. The key conclusions are:

\begin{itemize}
    \item the imaging of simulated ALMA Measurement Set with \texttt{resolve} results in a  successful recovery of the flux of the different components in the simulated dataset.
    \item In one of the first attempts to apply \texttt{resolve} on real ALMA data we obtain a high-fidelity image of protoplanetary disk Sz 114, highlighting the potential of \texttt{resolve} to create super-resolution images.
    \item For both test cases we obtain a robust estimation of the uncertainties on the measured fluxes, which is one of the major advantages of \texttt{resolve} compared with \texttt{tclean}.
\end{itemize}

With those encouraging results we highlight further areas where exciting developments can be made: \texttt{resolve} can be used to create spectral cubes, using correlation between channels as additional prior information; create combined maps of different antenna configuration, a particularly challenging case for \texttt{tclean} algorithm and where accurate assessment of the confidence in obtained data is especially necessary.

{\bf Acknowledgments:} This work made use of the following software: resolve \cite{Junklewitz.Bell.ea2015}, matplotlib \cite{Hunter2007}, NIFTy v.8 \cite{nifty}, astropy \cite{AstropyCollaboration.PriceWhelan.ea2018}, CASA \cite{McMullin.Waters.ea2007}, This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This paper makes use of the following ALMA data: ADS/JAO.ALMA\#2016.1.00484.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

\reftitle{References}


\section{Introduction} \label{introduction}

The chemically peculiar (CP) stars with
spectral types from early B to early F are traditionally characterized by the
presence of certain stellar absorption lines of abnormal strength or weakness that indicate
peculiar surface abundances \citep{1974ARA&A..12..257P}.
Current theories ascribe the observed chemical peculiarities to the interplay
between gravitational settling (atomic diffusion) and
radiative levitation \citep{2000ApJ...529..338R}.
Furthermore, convection, He settling, mass loss, and turbulent mixing can significantly
affect the chemical composition of CP stars \citep{2005A&A...443..627T}.

Many CP stars possess stable and globally organized magnetic fields with strengths of
up to several tens of kG \citep{2008AstBu..63..139R}. These stars are often referred
to as magnetic chemically peculiar (mCP) stars in the literature.
The origin of the magnetic field is still a matter of some
controversy \citep{2004Natur.431..819B}.
Recently, there was also classical pulsation among mCP stars found
\citep{2021mobs.confE..27H}.
Therefore, these objects are perfect astrophysical laboratories for investigating the most
important phenomena occurring on the upper main sequence in the presence of a stellar
magnetic field.

\citet{2020A&A...640A..40H} made use of the 520\,nm flux depression, typically for
mCP stars, to search for these objects by classifying spectra from the
Large Sky Area Multi-Object Fiber Spectroscopic
Telescope \citep[LAMOST,][]{2012RAA....12.1197C}. It
was shown \citep{2007A&A...469.1083K} that Fe is the principal contributor to
the 520\,nm depression for the whole range of  effective temperatures
($T_{\rm eff}$) of mCP stars,
while Cr and Si play a role primarily in the low $T_{\rm eff}$ region.
However, a strong 520\,nm flux depression was not found for all mCP stars
\citep{2006MNRAS.372.1804K}. The reason for this is most probably the
individual magnetic field
configuration together with the line of sight.

Motivated by the success of \citet{2020A&A...640A..40H}, we employed the $\Delta$a
photometric system which measures the above-mentioned flux depression
\citep{2005A&A...441..631P,2014A&A...562A..65S} to the low-resolution
blue photometer and red photometer (BP/RP) spectra of the $Gaia$ mission \citep{2021A&A...652A..86C}.
These spectra cover the wavelength region from 330 to 1050\,nm
with a resolving power between 25 and 100, depending on the
wavelength. The signal-to-noise ratio
 depends on the apparent magnitude and color of the object.
The Gaia consortium put a lot of effort into bringing all spectra
onto a common flux and pixel (pseudo-wavelength) scale, taking
into account variations over the focal plane. Finally, they
produced a mean spectrum from
all the observations of the same source \citep{2022arXiv220606205M}. The
spectra are available in a continuous representation, with a subset of
sources having spectra in a sampled form as well.
The $Gaia$ DR3 already includes about 219\,000\,000 mean BP/RP spectra
for objects up to the 18th magnitude in $G$.

The spectra are therefore perfectly suited to performing a case study if the known mCP stars
can be efficiently detected using the 520\,nm flux depression and synthetic
$\Delta$a magnitudes.

  \begin{figure}
   \centering
   \includegraphics[width=0.45\textwidth]{CMD.png}
      \caption{Diagram of  $M(G)_0$ vs $(BP-RP)_0$   for the mCP (upper panel) and
      normal-type stars (lower panel) together with isochrones from
      \citet{2012MNRAS.427..127B} for solar metallicity [Z]\,=\,0.0152\,dex. The two samples
        overlap.}
         \label{plot_CMD}
   \end{figure}

\section{Target selection} \label{target_selection}

For the apparent normal stars, we  chose all the objects included in the
papers by \citet{2005A&A...444..941P,2006A&A...458..293P}, who presented an
empirical $T_{\rm eff}$ calibration for B- to F-type stars. In total,
387 stars also have an available BP/RP spectrum. We note that
this sample also includes binary stars. Among them, 182 objects have a
``Dup'' flag equal to one and/or a RUWE parameter larger than 1.5, which are both
good indicators of binarity \citep{2022AJ....163...33Z}.

The mCP stars were chosen from the following references; some stars were listed
more than once:

\begin{itemize}
    \item \citet{2005A&A...441..631P}: We
     selected all stars with $\Delta$a\,$>$\,+20\,mmag from this catalogue of $\Delta$a measurements.
    \item \citet{2019ApJ...873L...5C}: They published a list of 157 mCP stars with
    resolved magnetically split lines from the Sloan Digital Sky Survey (SDSS)/Apache
    Point Observatory Galactic Evolution Experiment (APOGEE) survey.
    \item \citet{2020A&A...640A..40H}: Spectra from the LAMOST Survey
    showing a 520\,nm flux depression were used to find a sample of 1002 mCP stars.
\end{itemize}

From this large selection, 1240 mCP stars have a BP/RP spectrum. For both samples,
we made no cut according to the signal-to-noise ratios of the spectra.
To check if both samples are overlapping in the color--magnitude diagram (CMD),
we used the photometry and distances from the $Gaia$ DR3
\citep{2021AJ....161..147B,2022arXiv220800211G}. Furthermore, the reddening values
of the individual stars were taken from
\citet{2005A&A...444..941P,2006A&A...458..293P,2008A&A...491..545N,2020A&A...640A..40H}.
If no value was found, the reddening map of \citet{2019ApJ...887...93G} was applied.
The individual reddening values were transformed to the $Gaia$ filters according
to the recommendations by the $Gaia$
consortium\footnote{\url{https://www.cosmos.esa.int/web/gaia/edr3-extinction-law}}.
We note that most stars are in the solar neighborhood  ($d$\,$<$\,500\,pc),
and therefore the reddening is not significant.

In Fig. \ref{plot_CMD} the $M(G)_0$ versus $(BP-RP)_0$ diagram
is presented for the two samples,
together with isochrones from \citet{2012MNRAS.427..127B} for solar
metallicity [Z]\,=\,0.0152\,dex. The normal-type stars (including binaries)  cover the
whole spectral range and the main sequence very well. No stars can be found below the zero age
main sequence, which provides confidence in our method and reddening selection.
The mCP stars, except for a few outliers, are located in the same region in the CMD (i.e.,
in the same astrophysical parameter space).

     \begin{figure*}
   \centering
   \includegraphics[width=0.95\textwidth]{Delta_a.png}
      \caption{The $a$ versus ($g_1-y$) diagram for the 387 normal type (filled symbols)
      and 1240 mCP
          stars (open symbols). The normality line $a = -2.16(8) + 0.363(8) (g_1-y)$
          is defined as in the classical $\Delta$a photometric system. The dotted lines are
          the 95\% prediction bands to select mCP stars. We also indicated the approximate
          spectral types according to the CMD (Figure. \ref{plot_CMD}).}
         \label{plot_Delta_a}
   \end{figure*}

\begin{table}
\begin{center}
\caption{Statistics of the detection level for the synthetic $\Delta$a values
(Figure \ref{plot_Delta_a}). The numbers for ``+'' and ``$-$'' mean above and
below the 95\% prediction band of the normality line. The listed ratios
should   ideally be 100 for the
mCP stars and
0 for the normal-type stars.}
\label{result_classification}
\tiny
\begin{tabular}{cccc|ccc}
\hline
($g_1-y$) & $N_{\mathrm CP}-$ & $N_{\mathrm CP}$+ & \% &
$N_{\mathrm norm}-$ & $N_{\mathrm norm}$+ & \% \\
\hline
8.0$-$8.5 & 0 & 4 & 100 & 2 & 0 & 0 \\
8.5$-$9.0 & 3 & 25 & 89 & 32 & 0 & 0 \\
9.0$-$9.5 & 7 & 132 & 95 & 106 & 0 & 0 \\
9.5$-$10.0 & 12 & 202 & 94 & 85 & 0 & 0 \\
10.0$-$10.5 & 16 & 240 & 94 & 48 & 0 & 0 \\
10.5$-$11.0 & 15 & 247 & 94 & 25 & 0 & 0 \\
11.0$-$11.5 & 23 & 150 & 87 & 29 & 1 & 3 \\
11.5$-$12.0 & 13 & 82 & 86 & 26 & 0 & 0 \\
12.0$-$12.5 & 9 & 36 & 80 & 16 & 1 & 6 \\
12.5$-$13.0 & 10 & 14 & 58 & 7 & 0 & 0 \\
\hline
\end{tabular}
\end{center}
\end{table}

\section{Analysis} \label{analysis}

To calculate the synthetic magnitudes in the $\Delta$a photometric system, we
used the approach  published in \citet{2014A&A...562A..65S}.
We first normalized all spectra to the flux given at 402\,nm.
Then they were interpolated in the wavelength region from 480 to 580\,nm
to a one-pixel resolution of 0.1\,nm applying a standard polynomial technique.
This works perfectly because this wavelength region is very smooth.

For the filter curves,
the central wavelengths as published by \citet{1980A&A....83..328M}  were employed: $g_1$ (501\,nm),
$g_2$ (521.5\,nm), and $y$ (548.5\,nm) with a bandwidth of 13\,nm.
According to \citet{2014A&A...562A..65S}, this results in the most significant
$\Delta$a values.
The final magnitudes for the three filters are in arbitrary units. It is not
straightforward to determine the errors of the magnitudes. Because we did not limit
the signal-to-noise ratios (S/N) of the spectra, the source of this error was taken into
account. For each spectrum, we took the given S/N and
randomly added an error to the individual fluxes. For the lowest S/N, values of about 10, we
deduced an error of 1.5\%, which does not affect the conclusions and detection limits.

As the last step, we generated the $a$ versus ($g_1-y$) diagram for the normal type and mCP
stars, which is shown in Figure \ref{plot_Delta_a}. As can be seen, the normal stars
show a distinct band and are clearly separated from the mCP objects.
We determined the normality line as in the classical $\Delta$a photometric system (i.e.,
assuming that all stars exhibit the same interstellar reddening; peculiar objects deviate
from it by more than 3$\sigma)$. For the normal-type star sample, it was calculated as
\begin{equation}
    a = -2.16(8) + 0.363(8) (g_1-y)
\end{equation}
in the units of our synthetic photometric system. The characteristics are
similar to the observed normality
lines for open clusters
\citep{2007A&A...462..591N,2014A&A...564A..42P}. The color ($g_1-y$) is
well correlated with the
$T_{\rm eff}$, especially from late B- to early F-type stars.
We split the sample into ($g_1-y$) bins and counted how many stars are
below (``$-$'') and
above
(``+'') the normality line. The latter is the region where the mCP stars are
located. In Table \ref{result_classification} the results are presented for both
samples. In total, from 387 normal stars, only 2 are located slightly above the
prediction band of the region of the mCP objects.
From the 1240 mCP stars, 108 (8.7\%) are below the upper 95\% prediction band.
However, the distribution
of the latter is homogeneous over the ($g_1-y$) range. Most of them
are located among cooler-type objects. This effect is already known, as
\citet{2006A&A...448.1153K,2007A&A...469.1083K} showed that the capability
of the $\Delta$a photometric system drops significantly for stars with $T_{\rm eff}$
lower than 8000\,K. This is due to the observational fact that the 520\,nm
flux depression becomes less prominent compared to the overall line blanketing
\citep{2004MNRAS.352..863K}. For the B- and A-type stars we achieved detection
rates of almost 95\%, which is equavalent to a photometric accuracy of 10\,mmag
for classical CCD observations \citep{2005A&A...441..631P}.

\section{Conclusions} \label{conclusions}

We  showed that, with the help of the low-resolution BP/RP spectra of the
$Gaia$ mission, it is possible to efficiently search for and detect classical magnetic
chemically peculiar  (mCP) stars of the upper main sequence. This was done by
calculating magnitudes within the $\Delta$a photometric system which traces
the flux depression at 520\,nm typically for this star group. The latter are
benchmark objects for studying the effects of rotation, diffusion, mass loss,
and pulsation in the presence of an organized local stellar magnetic field.

Using well-known mCP stars and a sample of normal type stars (including about 50\% binaries),
an almost 95\% detection rate for B- and A-type stars was achieved. With our
technique it is therefore possible to pre-select the best candidates for further
spectroscopic and photometric studies.
Querying the different catalogues of the \textit{Gaia} DR3, we found spectral
tags indicating approximately 13 million B- and A-type stars up to $G$\,=\,17.65\,mag
among the published data release. \citet{2005MNRAS.362.1025P} estimated the lower
incidence of all CP stars in the Milky Way of 6\% with a published maximum value
of at least 16\%, respectively. The mCP stars account for about one-third of all
CP objects. Therefore, we adopt a mean value of 3\% for the
considered spectral type range. The ratio of mCP stars that are part of binary systems
is still inconclusive \citep{2020CoSka..50..570P}. In total, we estimate
$\sim 2 \times 10^{5}$ mCP stars to be potentially detectable for the discussed
accuracy of the presented synthetic $\Delta$a photometry.

\begin{acknowledgements}
This work was supported by the European Regional Development Fund,
project No. ITMS2014+: 313011W085.
This work presents results from the European Space Agency (ESA) space mission
Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis
Consortium (DPAC). Funding for the DPAC is provided by national institutions,
in particular the institutions participating in the Gaia MultiLateral Agreement
(MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia.
The Gaia archive website is https://archives.esac.esa.int/gaia.
This research has made use of the SIMBAD database,
operated at CDS, Strasbourg, France and
of NASA's Astrophysics Data System Bibliographic Services.
\end{acknowledgements}

\bibliographystyle{aa}


\section{Introduction}
Goal Recognition is the task of recognising the goal that an agent is trying to achieve from observations about the agent's behaviour in the environment \cite{DBLP:journals/firai/Van-HorenbekeP21,2018_geff_modelfree_modelbased}.
Typically, such observations consist of a trace (sequence) of executed actions in an agent's plan to achieve the goal, or a trace of world states generated by the agent's actions, while an agent goal is specified by a set of propositions.
Goal recognition has been studied in AI for many years, and it is an important task in several fields, including human-computer interactions \cite{human_interaction}, computer games \cite{min}, network security \cite{mirsky2019new}, smart homes \cite{harman2019action}, financial applications \cite{borrajo}, and others.

In the literature, several systems to solve goal recognition problems have been proposed \cite{DBLP:conf/ijcai/MeneguzziP21}.
The state-of-the-art approach is based on transforming a plan recognition problem into one or more plan generation problems solved by classical planning algorithms \cite{ramirez2009plan,pereira2020landmark,gr_noisy}.
In order to perform planning, this approach requires domain knowledge consisting of a model of each agent's action specified as a set of preconditions and effects, and a description of an initial state of the world, in which the agent perform the actions.
The computational efficiency (runtime) largely depends on the planning algorithm performance, which could be inadequate in a context demanding fast goal recognition (e.g., in real-time/online applications).\footnote{
Deciding plan existence in classical planning is PSPACE-complete  \cite{BYLANDER1994}.}

In this paper, we investigate an alternative approach
in which the goal recognition problem is formulated as a classification task, addressed through machine learning, where each candidate goal (a set of propositions) of the problem can be seen as a value class.
The primary aim is making goal recognition more accurate as well as faster by learning how to solve it in a given domain.
Given a planning domain specified by a set of propositions and a set of action names, we tackle the goal classification instances in the domain through a Recurrent Neural Network (RNN). A run of our RNN processes a trace of observed actions to compute how likely it is that each domain proposition is part of the agent's goal, for the problem instance under considerations. These predictions are then aggregated through a goal selection mechanism   to choose one of the candidate goals.

The proposed approach, that we call {\sf \small GRNet}, is generally faster than the model-based approach to goal recognition based on planning, since running a trained neural network can be much faster than plan generation. Moreover, {\sf \small GRNet}
 operates with minimal information, since the only information required as input for the trained RNN is a trace of action labels (each one indicating just the name of an observed action), and the initial state can be completely unknown.

The RNN is trained only once for a given domain, i.e., the same trained network can be used to solve a large set of goal recognition instances in the domain. On the other hand, as usual in supervised learning,  a  (possibly large) dataset of solved goal recognition instances for the domain under consideration is needed for the training.
When such data are unavailable or scarse, they can be synthetized via planning.
In such a case, the resulting overall system can be seen as a combined approach (model-based for generating training data, and model-free for the goal classification task) that outperforms the pure model-based approach in terms of both classification accuracy and classification runtime.
Indeed, an experimental analysis presented in the paper confirms that \ensuremath{{\sf GRNet}}\xspace achieves good performance in terms of both goal classification accuracy and runtime, obtaining consistently better performance with respect to a state-of-the-art goal recognition system over a class of benchmarks in six planning domains.

In rest of paper, after giving background and preliminaries about goal recognition and LSTM networks, we describe the \ensuremath{{\sf GRNet}}\xspace approach; then we present the experimental results; finally, we discuss related work and give the conclusions.

\section{Preliminaries}

We describe the problem goal of recognition, starting with its relation to activity/plan recognition, and we give the essential background on Long Short-Term Memory networks and the attention mechanism.

\label{sec:preliminaries}

\subsection{Activity, Goal and Plan Recognition}

Activity, plan, and goal recognition are related tasks \cite{geib2007plan}. Since
in the literature sometime they are not clearly distinguished, we begin with an informal definition of them following \cite{DBLP:journals/firai/Van-HorenbekeP21}.

Activity recognition concerns analyzing temporal sequences of (typically low-level) data generated by humans, or other autonomous agents acting in an environment, to identify the corresponding activity that they are performing.
For instance, data can be collected from wearable sensors, accelerometers, or images to recognize human activities such as running, cooking, driving, etc. \cite{DBLP:journals/firai/VrigkasNK15, JOBANPUTRA2019698}.

Goal recognition (GR) can be defined as the problem of identifying the intention (goal) of an agent from observations about the agent behaviour in an environment. These observations can be represented as an ordered sequence of discrete actions (each one possibly identified by activity recognition), while the agent's goal can be expressed either as a set of propositions or a probability distribution over alternative sets of propositions (each oven forming a distinct candidate goal).

Finally,
plan recognition is more general than GR and concerns both recognising the goal of an agent and identifying the full ordered set of actions (plan) that have been, or will be, performed by the agent in order to reach that goal; as GR, typically plan recognition takes as input a set of observed actions performed by the agent \cite{DBLP:journals/umuai/Carberry01}.

\subsection{Model-based and Model-free Goal Recognition}

In the approach to GR known as ``{goal recognition over a domain theory}'' \cite{geffner,DBLP:journals/firai/Van-HorenbekeP21,santos,gr_noisy},
the available knowledge consists of an underlying model of the behaviour of the agent and its environment. Such a model
represents the agent/environment states and the set of actions $A$ that the agent can perform; typically it is specified by a planning language such as
\textsc{pddl}\xspace \cite{mcdermott1998pddl}.
The states of the agent and environment are formalised as subsets of a set of propositions $F$, called {\em fluents} or {\em facts}, and each domain action in $A$ is modeled by a set of preconditions and a set of effects, both over $F$.  An instance of the GR problem in a given domain is then specified by:
\begin{itemize}
\item
an initial state $I$ of the agent and environment ($I \subseteq F$);
\item a sequence $O = \langle \textit{obs}_1, .., \textit{obs}_n \rangle$  of observations ($n \geq 1$), where each $obs_i$ is an action in $A$ performed by the agent;
\item  and a set $\mathcal{G} = \{ G_1,  .., G_m \}$ ($m \geq 1$) of possible goals of the agent, where each $G_i$ is a set of fluents over $F$.
\end{itemize}
The observations form a trace of the full sequence $\pi$ of actions performed by the agent to achieve the goal.
Such a plan trace is a selection of (possibly non-consecutive)  actions in $\pi$, ordered as in $\pi$.
Solving a GR instance consists of identifying the $G^*$ in $\cal G$ that is the (hidden) goal of the agent.

The approach based on a model of the agent's actions and of the agent/environment states, that we call {\em model-based goal recognition} (MBGR), defines GR as a reasoning task addressable by automated planning techniques \cite{DBLP:conf/ijcai/MeneguzziP21,GNT2016}.

An alternative approach to MBGR is {\em model-free goal recognition} (MFGR)  \cite{2018_geff_modelfree_modelbased,borrajo}. In this approach, GR
is formulated as a classification task addressed through machine learning.
The domain specification consists of a fluent set $F$, and a set of possible actions $A$, where each action $a\in A$ is specified by just a label (a unique identifier for each action).

A MFGR instance for a domain is specified by an observation sequence $O$
formed by action labels and, as in MBGR, a goal set $\cal G$ formed by subsets of $F$.
MFGR requires minimal information about the domain actions, and can operate without the specification of an initial state, that can be completely unknown.
Moreover, since running a learned classification model is usually fast, a MFGR system is expected to run faster than a MBGR system based on planning algorithms.
On the other hand, MFGR needs a data set of solved GR instances from which learning a classification model for the new GR instances of the domain.

\begin{example}

As a running example, we will use a very simple GR instance in the well-known \textsc{Blocksworld}\xspace domain. In this domain the agent has the goal of building one or more stacks of blocks, and only one block may be moved at a time. The agent can perform four types of actions:
{\tt \small Pick-Up} a block from the table, {\tt \small Put-Down} a block on the table,
 {\tt \small Stack} a block on top of another one, and {\tt \small Unstack} a block that is on another one.
 %
We assume that a GR instance in the domain involves at most 22 blocks. In \textsc{Blocksworld}\xspace
there are three types of facts (predicates):
{\small \tt On}, that has two blocks as arguments, plus {\small \tt On-Table} and {\small \tt Clear} that have one argument. Therefore, the fluent set $F$ consists of $22\times 21 + 22 + 22 = 506$ propositions.
The goal set ${\cal G}$ of the instance example consists of the two goals
$G_1=$ \{{\tt \small(On Block\_F Block\_C)}, {\tt \small (On Block\_C Block\_B)}\} and  $G_2=$ \{{\tt \small (On Block\_G Block\_H)}, {\tt \small (On Block\_H Block\_F)}\}; the observation sequence $O$ is $\langle${\tt \small (Pick-Up Block\_C)}, {\tt \small (Stack Block\_C Block\_B)}, {\tt \small (Pick-Up Block\_F)}$\rangle$.
\end{example}
\begin{figure*}
    \centering
    \includegraphics[scale=0.28]{imgs/architecture.png}
    \caption{
    Architecture of \ensuremath{{\sf GRNet}}\xspace. The input observations are encoded by embedding vectors and then fed to a LSTM neural network. After that the attention mechanism computes the context vector, which is used by a feed-forward layer to define the corresponding output values. This layer is composed by $|F|$ neurons, each one representing a possible fluent in the domain.
    The output of the neural network is then used by the instance component for selecting the goal with the highest score ($G_1$ in the example of the figure).
    The observed actions $a_{05}$, $a_{17}$, ..., $a_{22}$ are ordered from top to bottom according to their execution order.
    }
    \label{fig:architecture}
\end{figure*}

\subsection{LSTM and Attention Mechanism}
\label{sec:lstm}
A Long Short-Term Memory Network (LSTM) is a particular kind of Recurrent Neural Network (RNN). This kind of deep learning architecture is particularly suitable for processing sequential data like signals or text documents \cite{lstm}. With respect to the standard RNN, LSTM deals with typical issues such as vanishing gradient and long-term dependencies, obtaining better predictive performance \cite{Gers2000LearningTF}. Let $x_1, x_2 \ldots x_m$ be an input time series where $x_t \in \mathbb{R}^d$ is the feature vector representing the $t$-th element of the series, and $d$ is the dimension of each feature vector of the sequence. The long and short term memory states $c_{t}\in \mathbb{R}^N$  and $h_{t}\in \mathbb{R}^N$  at time step $t$ of the series, respectively, are computed recursively considering the values at the previous time step $t-1$  as follows:

\hspace{-0.2cm}
{\small
\begin{tabular}{ll}
$\hat{c}_t = \tanh(W_c [h_{t-1}, x_t] + b_c)$ &
$i_t = \sigma(W_i [h_{t-1}, x_t]  + b_i)$ \\
$f_t = \sigma(W_f [h_{t-1}, x_t]  + b_f)$&
$c_t = i_t * \hat{c}_t + f_t * c_{t-1} $\\
$o_t = \sigma(W_o [h_{t-1}, x_t] + b_o)$&
$h_t = \tanh(c_t) * o_t$\\
\end{tabular}
}
\noindent
where $\sigma$ denotes the sigmoid activation function and $*$ corresponds to the element-wise product; $W_f$, $W_i$, $W_o$, $W_c$ $\in \mathbb{R}^{(N+d)\times{N}}$ are the weight matrices and  $b_f$, $b_i$, $b_o$, $b_c$ $\in \mathbb{R}^N$ are the bias vectors; the vectors in square brackets are concatenated.
Weight matrices and bias vectors are typically initialized with the Glorot uniform initializer \cite{glorot2010understanding}, and they are shared by all the cells in the LSTM layer. $h_0$ and $c_0$ are initialized as zero vectors.

The {\it attention mechanism} \cite{DBLP:journals/corr/BahdanauCB14} is another layer which computes weights representing the contribution of each element of the sequence, and provides a
representation of the sequence (also called the \textit{context vector}) as the weighted average of the outputs ($h_t$) of the LSTM cells,
improving the predictive performance with respect to the base LSTM networks.
In our system, we use the so-called \textit{word attention} introduced by \citet{yang_attention} in the context of text classification.

\section{Goal Recognition through GRNet}

In this section we present our approach to goal recognition based on deep learning, \ensuremath{{\sf GRNet}}\xspace.
\ensuremath{{\sf GRNet}}\xspace is depicted in Figure \ref{fig:architecture}  consists of two main components.
The first component takes as input the observations of the GR instance to solve, and gives as output a score (between 0 and 1) for each proposition in the domain proposition set $F$. This component,  called {\it Domain Component}, is general in the sense that it can be used for every GR instance over $F$ (training is performed once for each domain).
The second component, called {\em Instance Component}, takes as input the proposition ranks generated by the domain component for a GR instance, and uses them to select a goal from the candidate goal set ${\cal G}$.

\subsection{The Domain Component of GRNet}
Given a sequence of observations, represented on  the left side of Figure \ref{fig:architecture}, each action $a_i$ corresponding to an observation is encoded as a vector $e_i$ of real numbers by an embedding layer \cite{bengio2003neural}.\footnote{\url{https://keras.io/api/layers/core_layers/embedding/}}
In Figure \ref{fig:architecture}, the observed actions are displayed from top to bottom in the order in which they are executed by the agent. The embedding layer is initialised with random weights, and trained at the same time with the rest of the domain component.

The index of each observed action is simply the result of an arbitrary order of the actions that is computed in the pre-processing phase, only once for the domain under consideration.
Please note that two actions $a_i$ and $a_j$ consecutively observed may not be consecutive actions in the full plan of the agent
(the full plan may contain any number of actions between $a_i$ and $a_j$).

The sequence of embedding vectors is then fed to a LSTM Neural Network, and the result of the output of each cell is processed by the Attention Mechanism.
After computing a weight for the contribution of each cell, this layer provides a so-called \textit{context-vector}
that summarizes the information contained in the trace of plan.
The context vector is then passed to a feed-forward layer, which has $N$ output neurons with \textit{sigmoid} activation function.
$N$ is the number of the domain fluents (propositions) that can appear in any goal of {$\cal G$} for any GR instance in the domain; for our experiments $N$ was set to the size of the domain fluent set $F$, i.e., $ N = |F|$.
The output of the $i$-th neuron $\overline{o}_i$ corresponds to the $i$-th fluent $f_i$ (fluents are lexically ordered), and the activation value of $\overline{o}_i$ gives a rank for $f_i$ being true in the agent's goal (with rank equal to one meaning that $f_i$ is true in the goal).
In other words, our network is trained as a multi-label classification problem, where each domain fluent can be considered as a different binary class. As loss function, we used standard binary crossentropy.

As shown in Figure \ref{fig:architecture}, the dimension of the input and output of our neural networks depend only on the selected domain and some basic information, such as the maximum number of possible output facts
that we want to consider.
The dimension of the embedding vectors, the dimension of the LSTM layer and other hyperaparameters of the networks are selected using the Bayesian-optimisation approach provided by the Optuna framework \cite{akiba2019optuna}, with a validation set formed by 20\% of the training set while the remaining 80\% is used for training the network.
More details about the hyperparameters are given in the Supplementary Material.

\subsection{The Instance Specific Component of GRNet}
After the training and optimisation phases of the domain component, the resulting network can be used
 to solve any goal recognition instance in the domain through the instance-specific component of our system (right part of Figure \ref{fig:architecture}).
Such component performs an evaluation of the candidate goals in $\cal G$ of the GR instance, using the output of the domain component fed by the observations of the GR instance.
To choose the most probable goal in $\mathcal{G}$ (solving the multi-class classification task associated to the GR instance), we designed a simple score function that indicates how likely it is that $G$ is the correct goal, according to the neural network of the domain component.
This score is defined as
$    score(G, \overline{o}) = \sum_{f \in G} \overline{o}_f
$,
where $\overline{o}$ is the network output vector, and  $\overline{o}_f$ is the network output for fact $f$.
For each candidate goal $G \in \cal G$, we consider only the output neurons that have associated facts in $G$, ignoring the others. By summing only these predicted values we can derive an overall score for $G$ being the correct goal.
We compute this score for all candidate goals in $\cal G$, and we select the one with the highest score ($G_1$ in Figure \ref{fig:architecture}).

\begin{example}
In our running example, since we are assuming that GR instance involve at most 22 blocks, we have that the action set $A$ is formed by
$22$ {\small \tt Pick-Up} actions,  22 {\small \tt Put-Down} actions, $22 * 21 = 462$ {\small \tt Stack} actions and  462 {\small \tt Unstack} actions, for a total of  $968 = |A|$ different actions in the domain.
Suppose that the three observed actions
{\small \tt (Pick-Up Block\_C)}, {\small \tt (Stack Block\_C Block\_B)} and {\small \tt (Pick-Up Block\_F)}
forming the observation sequence $O$
have ids corresponding to indices $5$, $17$ and $21$, respectively.
In the Domain Component of GRNet, after being processed by the embedding layer, the input $O$ is represented by the sequence of vectors $e_{05}$, $e_{17}$ and $e_{21}$.
Then this sequence is fed to the LSTM layer and subsequently to the attention mechanism, producing a context vector $c$ representing the entire plan trace formed
by the observed actions. Finally, vector $c$ is
processed by a final feed-forward layer made of
$\abs{F} = 506$ neurons.
In this representation, each neuron corresponds to a distinct proposition in $F$. Therefore, if the network has to predict candidate goal $G_1$, made by
{\small \tt (On~Block\_C~Block\_B)} and {\small \tt (On~Block\,\,F~Block\_C)},
their corresponding neuron should have value 1, while the  neurons of the different propositions in $G_2$ should
have value zero.
Our GR instance is made by the plan trace {\small \tt (Pick-Up Block\_C)}, {\small \tt (Stack Block\_C Block\_B)} and {\small \tt (Pick-Up Block\_F)} and by two possible goals: $G_1$, made by {\small \tt (On Block\_C Block\_B)} and {\small \tt (On Block\_F Block\_C)}, and $G_2$ made by  {\small \tt  (On Block\_G Block\_H)} and {\small \tt (On Block\_H Block\_F)}.
Therefore, in the Instance Component of GRNet, we calculate the prediction values of $G_1$ and $G_2$ as the sum of the predictions for the neurons representing their facts.
Suppose that $\overline{o}_{{\tt (On~Block\_C~Block\_B)}}=1.000$, $\overline{o}_{{\tt  (On~Block\_F~Block\_C)}}=0.017$, $\overline{o}_{{\tt   (On~Block\_G~Block\_H)}}=0.000$, $\overline{o}_{{\tt (On~Block\_H~ Block\_F)}}=0.003$, we have that the final score of $G_1$ is $1.017$, while the final score of $G_2$ is $0.003$.
The goal with the highest score ($G_1$) is selected as the most probable goal solving the GR instance.
\end{example}

\section{Experimental Analysis}

After describing the benchmark domains, datasets, and GR instances, we
analyse the performance of \ensuremath{{\sf GRNet}}\xspace, comparing it with a state-of-the-art goal recognition system \cite{pereira2020landmark}.

\subsection{Benchmark Suite and Data Sets}
\label{sec:setup}
\begin{table}
\scriptsize
    \centering
    \begin{tabular}{ lccccc }

        \toprule
        Domain & $\abs{A}$ & $\abs{F}$ & $\abs{G_i}$ & $\abs{\mathcal{G}}$ \\
        \midrule
        \textsc{Blocksworld}\xspace  & 968 & 506  &   [4,16] & [19,21]\\
        \textsc{Depots}\xspace       & 13050 & 150 &   [2,8]  & [7,10]\\
        \textsc{Driverlog}\xspace    & 4860 & 156 &   [4,11] & [6,10]\\
        \textsc{Logistics}\xspace    & 15154 &  154 & [2,4]  & [10,12]\\
        \ensuremath{\textsc{Satellite}}\xspace    & 33225 & 629  & [4,9]  & [6,8]\\
        \textsc{Zenotravel}\xspace   & 23724 & 66  &  [5,9]  & [6,11]\\
        \bottomrule
    \end{tabular}
    \caption{Size of $A$, $F$, $G_i \in \cal G$ and $\cal G$ in the considered GR instances for each considered domain.  Interval $[x,y]$ indicates a range of integer values from $x$ to $y$.}
    \label{tab:experiments_detail}
\end{table}
 We  consider six benchmark domains that are well known in the automated planning community: \textsc{Blocksworld}\xspace, \textsc{Depots}\xspace, \textsc{Driverlog}\xspace, \textsc{Logistics}\xspace, \ensuremath{\textsc{Satellite}}\xspace and \textsc{Zenotravel}\xspace \cite{ipcfirst,ipcthird}.

In order to create the (solved) GR instances for the training and test sets in the considered domains, we used automated planning techniques.
Concerning the training set,
for each domain, we randomly generated a large collection of (solvable) plan generation problems
of different size.
We considered the same ranges of the numbers of involved objects as in the experiments of  \citet{pereira2020landmark}.
For each of these problems, we computed up to four (sub-optimal) plans solving them. As planner we used  \textsc{lpg}\xspace \cite{lpg_jair}, which allows to specify the number of requested different solutions for the planning problem it solves.
From the generated plans we derived the observation sequences for the training samples by randomly selecting actions from the plans (preserving their relative order). The selected actions are between
30\% and 70\% of the plan actions.
The generated training set consists of pairs $(O, G^*)$ where $O$ is a sequence of observed actions obtained by sampling a plan $\pi$, and $G^*$ is the hidden goal corresponding to the goal of the planning problem solved by $\pi$.
For each considered domain, we created a training set with $55000$ pairs.

For evaluating \ensuremath{{\sf GRNet}}\xspace we generated two test sets formed by GR instances {\em not seen at training time}: \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace and \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace.
Such test instances were generated as for the train instances, except that the observation sequences were derived using also plans generated by a planner different from the one used for creating the train instances (\textsc{lpg}\xspace). As second planner we used \textsc{lama}\xspace \cite{richter:10:lama}.

 \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace is a generalisation and extension of the test set used in \cite{pereira2020landmark} for the same domains that we consider. We indicate this original test set with \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace. \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace includes all
 \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace instances, and the goal sets ($\mathcal{G}$) of \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace and \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace are the same also for the new instances.
The additional instances in \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace are motivated by the limited number and particular structure of those in \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace.
The structure is limited because
the observations in the instances of \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace are created from plans for the goals in $\mathcal{G}$ that are all generated from the {\em same} initial state. In \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace, the GR instance are created  considering multiple plans generated from different initial states, obtaining a richer diversification of the observation traces and a larger number of test instances.
 In particular, for each of \textsc{Depots}\xspace, \textsc{Driverlog}\xspace, \ensuremath{\textsc{Satellite}}\xspace and \textsc{Zenotravel}\xspace, \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace contains only $84$ instances, while \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace test set contains roughly $6000$ instances.

For each plan generated for being sampled, we randomly derived three different action traces formed by 30\%,  50\% and 70\% of the plan actions, respectively. This gives three groups of test instances, for each considered domain, allowing to evaluate the performance of \ensuremath{{\sf GRNet}}\xspace also in terms of different amounts of available observations.
Table \ref{tab:experiments_detail}  gives information about the size of the GR instances in our test and training sets for each domain, in terms of number of possible actions ($\abs{A}$), facts ($\abs{F}$),  min/max size of a goal ($\abs{G_i}$) in a goal set $\cal G$, and min/max size of a goal set ($\abs{\mathcal{G}}$). More details in the supplementary material.

Test set \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace was created to evaluate
how well the compared systems behave on GR instances of different difficulty. We focus this analysis on a specific domain (\textsc{Zenotravel}\xspace).
In \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace, the generated GR instances are grouped into several classes according to their difficulty. As difficulty measure
we used the notion of \emph{recognizability of the hidden goal}, which is inspired by the notion of the ``uniqueness of landmarks" introduced by \citet{pereira2020landmark}.
Specifically, the recognizability $R(G)$ of a goal $G \in \cal G$ is defined as

$R(G) = \sum_{f \in G}{\frac{1}{|\{G' \; \mid \; G' \in \mathcal{G} \land f \in G'\}|}}$.

\noindent
The lower $R(G)$ is, the more difficult recognising $G$ is; vice versa, the higher $R(G)$ is, the more discernible $G$ is.

We normalized $R_Z(G)$ as a value between $0$ and $1$, denoted $R_Z(G)$.
For example, if $\mathcal{G} = \set{G_1, G_2, G_3 }$, with $G^* = G_1 = \set{a,b,c}$, $G_2 = \set{ a,e,f }$ and $G_3 = \set{g,h,i}$, then $R(G^*)= \frac{1}{2} + 1 +1 = \frac{5}{2}$ and $R_Z(G^*) = 0.75$ (high recognizability). While, if
$\mathcal{G} = \set{G_1, G_2, G_3}$ with $G^* =  G_1 = \set{a,b,c}$, $G_2 = \set{a, b, x}$ and $G_3 = \set{a,b,y}$, then $R(G^*)=\frac{1}{3} + \frac{1}{3} + 1 = \frac{5}{3}$, and so $R_Z(G^*) = 0.33$ (low recognizability).

Using different values for $R_Z(G^*)$, we have generated nine classes of GR instances, denoted $C_1,...,C_9$. For each GR instance in class $C_i$, we have $0.1 \cdot i \leq R_Z(G^*) < 0.1 \cdot (i+1)$, for $i = 1...9$. Each class consists of 100 GR instances.

\subsection{Experimental Results}
\label{sec:results}

We present the results obtained by \ensuremath{{\sf GRNet}}\xspace and the state-of-the-art system {\ensuremath{\sf LGR}}\xspace by \citet{pereira2020landmark}  over the test sets \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace and \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace. In order to have a fair comparison, for {\ensuremath{\sf LGR}}\xspace we used the $h_{\textit{uniq}}$  heuristic for the goal selection, because the authors stated that it performs better than others. In {\ensuremath{\sf LGR}}\xspace it is possible to set a threshold $0\leq \theta \leq 1$ that controls the confidence to accept a goal as a possible solution for the given GR instance. In order to select a single goal in $\cal G$, we set $\theta=0$. Nevertheless, we observed that in some GR tasks {\ensuremath{\sf LGR}}\xspace returns more than one goal. %
As done in \cite{pereira2020landmark}, for {\ensuremath{\sf LGR}}\xspace  we considered solved a GR instance when their system returns a set of goals {\em containing} the correct goal (making the evaluation more favourable for {\ensuremath{\sf LGR}}\xspace).

Goal recognition {\em accuracy} for a set of test instances is defined as the percentage of instances whose goals are correctly identified (predicted) over the total number of the instances in the set.

\begin{table}[]
    \centering
    \scriptsize
    \begin{tabular}{ l|cc|cc|cc }
        \toprule
        \multirow{2}{*}{Domain} & \multicolumn{2}{c |}{30\% of  the plan} & \multicolumn{2}{c|}{50\% of the plan} & \multicolumn{2}{c}{70\% of the plan} \\
        & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace\\
        \midrule
        \textsc{Blocksworld}\xspace & 38.9 & \bf{53.4*} & 52.9 & \bf{69.4*} & 76.0 & \bf{82.6} \\
        \textsc{Depots}\xspace & 45.9 & \bf{60.9*} & 65.4 & \bf{75.8*} & 82.1 & \bf{86.3} \\
        \textsc{Driverlog}\xspace & 43.9 & \bf{61.4*} & 61.5 & \bf{74.8*} & 79.7 & \bf{83.9} \\
        \textsc{Logistics}\xspace & 50.1 & \bf{65.7*} & 67.8 & \bf{78.5*} & 82.6 & \bf{87.4} \\
        \ensuremath{\textsc{Satellite}}\xspace & 69.0 & \bf{71.9} & 83.2 & \bf{84.0} & 90.3 & \bf{93.7} \\
        \textsc{Zenotravel}\xspace & 47.9 & \bf{77.0*} & 68.8 & \bf{89.7*} & 87.1 & \bf{96.0} \\
        \bottomrule
    \end{tabular}
    \caption{Goal classification accuracy (percentages of GR instances correctly solved) of {\ensuremath{\sf LGR}}\xspace and \ensuremath{{\sf GRNet}}\xspace over six domains. Bold results indicate better performance;
    ``*" indicates accuracy improvements of \ensuremath{{\sf GRNet}}\xspace wrt {\ensuremath{\sf LGR}}\xspace of at least $10$ points.
    }
    \label{tab:results}
\end{table}

\paragraph{Results for \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace.}
Table \ref{tab:results} summarizes the performance results of {\ensuremath{\sf LGR}}\xspace and \ensuremath{{\sf GRNet}}\xspace in terms of accuracy with test set \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace.
We focus the analysis on test instances derived from a mix of plans generated by \textsc{lama}\xspace  and \textsc{lpg}\xspace (half of the instances from each of the planners, for every domain and percentage of observed actions). It should be noted that very similar performance results were obtained using plans from either only \textsc{lpg}\xspace or only \textsc{lama}\xspace; these results are in the supplementary material.

\ensuremath{{\sf GRNet}}\xspace performs generally well.
Even for instances with only $30\%$ of the actions, it shows interesting performance, with the best result for \textsc{Zenotravel}\xspace where \ensuremath{{\sf GRNet}}\xspace reaches accuracy $77$.
With 50\% of the actions, the accuracy of \ensuremath{{\sf GRNet}}\xspace improves in every domain by more than $10$ points. For instance, in $\ensuremath{\textsc{Satellite}}\xspace$ \ensuremath{{\sf GRNet}}\xspace improves from $72.8$ to $84.0$. With 70\% of the actions, our system improves further reaching accuracy above $80$ in all domains, and obtaining an impressive performance for \textsc{Zenotravel}\xspace ($95.9$).
Compared to {\ensuremath{\sf LGR}}\xspace, \ensuremath{{\sf GRNet}}\xspace performs always better, for every considered domain and percentage of observations.
In many cases, the performance improvement wrt {\ensuremath{\sf LGR}}\xspace is at least of $10$ points.

Moreover, \ensuremath{{\sf GRNet}}\xspace's performance does not seem to be affected by the diversity of the domains indicated by the four parameters of Table \ref{tab:experiments_detail}.
While the remarkable performance obtained for \textsc{Zenotravel}\xspace might be correlated with the fact that in this domain the test instances have only $66$ facts (see column $\abs{F}$ of Table \ref{tab:experiments_detail}),  the results for \ensuremath{\textsc{Satellite}}\xspace are not so distant even if the instances in this domain have $629$ facts. Analysing the experimental results, it seems that also the number of the actions has no significant impact on the performance. In fact, while \textsc{Blocksworld}\xspace has only $968$ actions, the other domains have more than $15000$ actions, and \ensuremath{{\sf GRNet}}\xspace obtain better results for them. This is probably due to the embedding layer that is able to learn a compact and informative representation even with a large vocabulary of actions. Overall, \ensuremath{{\sf GRNet}}\xspace exhibits good robustness with respect to the size of the space of actions and the number of facts in the domains (the output of the network).

In term of CPU time required to solve (classify) a GR instance, \ensuremath{{\sf GRNet}}\xspace performs generally much better than {\ensuremath{\sf LGR}}\xspace. The average execution time of {\ensuremath{\sf LGR}}\xspace is $1.158$ seconds with a standard deviation of $0.87$ seconds, while \ensuremath{{\sf GRNet}}\xspace runs on average in $0.06$ seconds with a standard deviation of $0.04$ seconds. For lack of space, details about CPU times for each domain are given in the Supplementary Material.

While we consider the evaluation of \ensuremath{{\sf GRNet}}\xspace using \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace more significant and more informative than using  \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace, we performed a comparison of
\ensuremath{{\sf GRNet}}\xspace and {\ensuremath{\sf LGR}}\xspace also using \ensuremath{\textsc{TS}_{\textit{Per}}}\xspace, obtaining results that overall are in favor of \ensuremath{{\sf GRNet}}\xspace also for this restricted test set.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.90\linewidth]{imgs/recogn.png}
    \caption{
    Accuracy results of {\ensuremath{\sf LGR}}\xspace and \ensuremath{{\sf GRNet}}\xspace  on several GR instances grouped into classes of decreasing difficulty with test set \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace. $C_1$ is the most difficult class while $C_9$ is the easiest one.}
    \label{fig:uniq_exp}
\end{figure}

\paragraph{Results for \ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace.}
Figure \ref{fig:uniq_exp} shows the accuracy of the two compared systems considering different classes of test sets with decreasing difficulty measured using $R_Z$. As expected, the accuracy of \ensuremath{{\sf GRNet}}\xspace depends on the difficulty of the problem, since there is an increasing trend in terms of accuracy for each observation percentage. This trend is more evident when we have $30\%$ of the actions and becomes less marked as the number of observations grows. {\ensuremath{\sf LGR}}\xspace appears to be more stable over the recognizability classes than \ensuremath{{\sf GRNet}}\xspace. However, \ensuremath{{\sf GRNet}}\xspace always performs significatly better than {\ensuremath{\sf LGR}}\xspace regardless the value of $R_Z$.

\begin{figure}
    \centering
    \includegraphics[width=0.64\linewidth]{imgs/incremental_our.png}
    \caption{Accuracy of \ensuremath{{\sf GRNet}}\xspace trained using data sets of GR examples with different sizes (percentages of the original train set) using test set \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace
    in domain \ensuremath{\textsc{Satellite}}\xspace.}
    \label{fig:incr_training_our}
\end{figure}

\paragraph{Sensitiveness to the training set size.}
Since the predictive performance of a machine learning system can be deeply influenced by the number of training instances, we experimentally investigated how much \ensuremath{{\sf GRNet}}\xspace is sensible to the this issue. We focused this analysis on domain \ensuremath{\textsc{Satellite}}\xspace, training several neural networks with different fractions of our training set: 20\%, 40\%, 60\% and 80\%. Figure \ref{fig:incr_training_our} shows how accuracy increases for
\ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace
when the training set size increases. In particular, for \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace we can observe that using only 20\% of the training instances gives accuracy lower than 40 in all three cases (30-50-70\% of observed actions), but accuracy rapidly improves reaching more than 60 using 60\% of the training instances.

We evaluated \ensuremath{{\sf GRNet}}\xspace also for larger training sets, up to twice the number of instances in the original training set.
As it can be seen in Figure \ref{fig:incr_training_our},
the enlarged training set for
\ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace produces only a small improvement in accuracy.

\section{Related Work}
Goal recognition and plan recognition have been extensively studied through model-based approaches exploiting planning techniques \cite{DBLP:conf/ijcai/MeneguzziP21,geffner,gr_noisy,santos,pereira2020landmark} or matching techniques relying on plan libraries (e.g., \cite{plan_libraries, mirsky}).
We have presented a model-free approach,
discussing its advantages and experimentally comparing it with the state-of-the-art model-based approach {\ensuremath{\sf LGR}}\xspace \cite{pereira2020landmark}.
This is a planning-based system exploiting landmarks \cite{hoffmann2004ordered}.
Differently from \ensuremath{{\sf GRNet}}\xspace, {\ensuremath{\sf LGR}}\xspace
 requires domain knowledge and performs no learning from previous experiences. Similarly to \ensuremath{{\sf GRNet}}\xspace , the output goal provided by {\ensuremath{\sf LGR}}\xspace is not guaranteed to be correct.

Many approaches to human activity recognition adopt architectures similar to the RNN of \ensuremath{{\sf GRNet}}\xspace \cite{9760794, yin2022novel, chen2016lstm, khatun2022deep}. These systems, though, address a problem different from goal recognition; they deal with noisy input data from sensors, and perform a specific classification task (with fixed classification values). As a consequence, these architectures provide solutions to very specific problems. The work in this paper addresses the goal recognition problem, it deals with lack of observability in the actions of the agent's plan, and proposes a more general approach allowing to solve, by the same trained network, different goal recognition instances in the domain.

Concerning GR systems using neural networks, some works use them
for specific applications, such as game playing \cite{min}. \ensuremath{{\sf GRNet}}\xspace is more general, as it can be applied to any domain of which we know sets $F$ and $A$.
In order to extract useful information from image-based domains and perform goal recognition, \citet{amado_images} used a pre-trained encoder and a LSTM network for representing and analysing sequence of observed states, rather than actions as in our approach.
\citet{amado_observations} trained a LSTM-based system to identify missing observations about states in order to derive a more complete sequence of
states by which a MBGR system can obtain better performance.

\citet{borrajo} investigated the use of XGBoost \cite{chen2016xgboost} and LSTM neural networks for goal recognition using
only traces of plans, without any knowledge on the domain, similarly to  our approach.
However, \citeauthor{borrajo} train a specific machine learning model for each goal recognition instance (the goal set $\cal G$ is fixed), using instance-specific datasets for training and testing. Instead, in our approach we train a general-purpose neural network that can be used to solve a large number of different goal recognition instances, without the need of designing or training a new model.
Moreover, the experimental evaluation of the networks proposed in \cite{borrajo}
use peculiar goal recognition benchmarks with custom-made instances.
Instead, in our work we evaluate \ensuremath{{\sf GRNet}}\xspace much more in depth using known benchmarks introduced by \cite{pereira2020landmark}, an extension of them having many more testing instances, and additional benchmarks based with different degrees of goal recognition difficulty.

\citet{maynard} compared model-based techniques and approaches based on deep learning for goal recognition. However, as in \cite{borrajo}, such a comparison is made using specific instances, and several kinds of neural networks are trained to
predict directly the goal among a set of possible ones, instead of the facts that belong to the goal as in our approach. This makes the trained networks in \cite{maynard} specific for the considered GR instances in a domain, while our approach is more general since it trains a single network for the  domain.

Another substantial difference is that, while in a typical goal recognition problem we can have missing observation across the entire plan of the agent(s), the work in \cite{maynard} considers only observations from the start of the plan to a given percentage of it, treating every possible successive observation as missing.

\citet{amado2022goal} proposed a framework that combines off-the-shelf model-free reinforcement learning and state-of-the-art goal recognition techniques achieving promising results. However, similarly to \cite{borrajo}, their approach is designed to solve a specific goal recognition instance where the goal set $\mathcal{G}$ is fixed. On the contrary, the trained RNN of \ensuremath{{\sf GRNet}}\xspace can  be used to solve all GR instances definable over the domain and action sets ($F$ and $A$).

\section{Conclusions}

We have proposed an approach to address goal recognition as a deep learning task. Our system, called \ensuremath{{\sf GRNet}}\xspace, learns to solve (classify) goal recognition tasks from past experience in a given domain, and
requires no model of the domain actions nor a specification of an initial state.
Learning consists in training only one neural network for the considered domain, allowing to solve a large collection of GR instances in the domain by same trained network.
An experimental analysis shows that \ensuremath{{\sf GRNet}}\xspace performs well in several benchmark domains,
in terms of both accuracy and runtime of the trained system, outperforming a state-of-the-art GR system.

The GR tasks addressable by \ensuremath{{\sf GRNet}}\xspace are limited to those involving subsets of  facts and actions appearing in the training set. An interesting question for future work is how to extend \ensuremath{{\sf GRNet}}\xspace to solve GR instances involving new actions and facts.
We also intend to (i) carry out more experiments to evaluate \ensuremath{{\sf GRNet}}\xspace and
investigate how its RNN can rich high performance (e.g., by analysing the weights provided by the attention mechanism and the learning process), and (ii) study the use of the more recent deep learning architectures.

Finally, an interesting direction for future work concerns finding effective ways of combining the model-free and the model-based approaches, the first based on learning from past experience and the second exploiting automated reasoning.
While a loose combination is using planning to generate samples in the training dataset for a learning-based GR system, finding a tighter integration that can further improve goal recognition performances is for future research.

\newpage

\section{Supplementary Material}
\label{sec:supplementary-data}

\subsection{Domains description}
We provide a very brief description of the six well-known planning domains used in our experiments:
\begin{itemize}
\item \textbf{\textsc{Blocksworld}\xspace}. The domain consists of an hand robot that has to stack or unstack blocks, picking up them one at a time, in order to obtain a desired configuration of an available set of blocks.

\item \textbf{\textsc{Depots}\xspace}. The domain consists of actions for loading and unloading packages into trucks through hoists, and moving them between depots. The goals concern having the packages at certain depots.

\item \textbf{\textsc{Driverlog}\xspace}. In this domain there are drivers (trucks) that can walk (drive) between locations. Walking and driving require traversing different paths. Packages can be loaded into or unloaded from trucks, that can be moved by drivers. The goal is to deliver (move) all packages (drivers) to their destinations.

\item \textbf{\textsc{Logistics}\xspace}. In this domain there are  aircrafts that can fly between cities, trucks that can move between locations within a city, and packages that can be loaded into/unloaded from trucks and aircrafts. The goal is to deliver a set of packages to their delivery locations.

\item \textbf{\ensuremath{\textsc{Satellite}}\xspace}. This is a scheduling domain in which one or more satellites can make certain observations, collect data, and download the collected data to a ground station. The goals concern having observation data at a ground station.

\item \textbf{\textsc{Zenotravel}\xspace}. This is another variant of a  transportation domain where passengers have to be embarked and disembarked into aircrafts, that can fly between cities at two possible speeds. The goals concern transporting (move) all passengers (aircrafts) to their required destinations.
\end{itemize}

\subsection{Size of the GR instances for the training/test sets in the considered benchmark domains}\label{ranges}
We give details about the number of objects involved in the GR instances.
For each object type in a domain, we report its name, the minimum and the maximum number of objects of that type ($min$ and $max$) that are involved in an instance, and the total number of objects of that type in the domain ($objs$); $objs$ indicates the number of all possible objects of a type that can be used to define a GR instance using at most $max$ of them. We chose these $min$-$max$ ranges because they are the {\em same} as those used in the experiments of \citet{pereira2020landmark}.

\begin{itemize}
    \item \textsc{Blocksworld}\xspace. \set{
    block: \set{min: 7, max: 17, objs: 22}
    }

    \item \textsc{Depots}\xspace, \set{
    depot: \set{min: 1, max: 3, objs: 3 },\;
    distributor: \set{min: 1, max: 3, objs: 3},\;
    truck: \set{min: 2, max: 3, objs: 3 },\;
    pallet: \set{min: 2, max: 6, objs: 6},\;
    crate: \set{min: 2, max: 10, objs: 10},\;
    hoist: \set{min:2, max: 6, objs: 6 }
    }

    \item \textsc{Driverlog}\xspace, \set{
    driver: \set{min: 2, max: 3, objs: 3},\;
    truck: \set{min: 2, max: 3, objs: 3},\;
    package: \set{min: 2, max: 7, objs: 7},\;
    locations_s: \set{min: 3, max: 12, objs: 12},\;
    location_p: \set{min: 2, max: 25, objs: 41}
    }

    \item \textsc{Logistics}\xspace, \set{
    airplane: \set{min: 1, max: 8, objs: 8},\;
    airport: \set{min: 2, max: 8, objs: 8},\;
    location: \set{min: 6, max: 11, objs: 11},\;
    city: \set{min: 2, max: 6, objs: 6},\;
    truck: \set{min: 2, max: 5, objs: 5},\;
    package: \set{min: 2, max: 14, objs: 14}
    }

    \item \ensuremath{\textsc{Satellite}}\xspace, \set{
    satellite: \set{min: 1, max: 5, objs: 5},\;
    instrument: \set{min: 1, max: 11, objs: 11},\;
    mode: \set{min: 3, max: 5, objs: 12},\;
    direction: \set{min: 7, max: 17, objs: 37}}

    \item \textsc{Zenotravel}\xspace, \set{
    aircraft: \set{min: 2, max: 3, objs: 3},\;
    person: \set{min: 5, max: 8, objs: 8},\;
    city: \set{min: 3, max: 6, objs: 6},\;
    flevel: \set{min: 7, max: 7, objs: 7}
    }

\end{itemize}

\subsection{Hyperparameters of the Neural Network}

Table \ref{tab:hyperparameters_range_sup} reports the hyperparameters range used in the Optuna framework. For each benchmark domain, we used a separate process of optimization (study) which executes 30 objective function evaluations (trials). We used a sampler that implements the Tree-structured Parzen Estimator algorithm.

\vspace{0.2cm}
\noindent
Table \ref{tab:hyperparameters_sup} reports the hyperparameters of the neural networks in our experiments. For all the experiments we selected 64 as batch size; we used Adam as optimizer with $\beta_1 = 0.9 $ and $\beta_2 = 0.99$.

\newpage

\begin{table}[h!]
\centering

\begin{tabular}{lc}
\toprule
Hyperparameter & Range\\
\midrule
${\abs{E}}$ & [50, 200]\\
${\abs{LSTM}}$ & [150, 512]\\
Use Dropout & \{True, False\}\\
Dropout & [0, 0.5]\\
Use Rec. Dropout & \{True, False\}\\
Rec. Dropout & [0, 0.5]\\

\bottomrule
\end{tabular}
\caption{
Value ranges of the hyperparameters for the Bayesian-optimisation done by the Optuna framework. ${\abs{E}}$ represents the dimension of the embedding vectors, ${\abs{LSTM}}$ is the number of neurons in the LSTM layer.  Interval [x, y] indicates a range of integer values from x to y, while set $\{x_1,..x_n\}$ enumerates all possible values the hyperparameter can assume.}
\label{tab:hyperparameters_range_sup}
\end{table}

\begin{table}[h!]
\centering

\begin{tabular}{lcccc}
\toprule
Domain & ${\abs{E}}$ & ${\abs{LSTM}}$ & Dropout & Recurrent Dropout\\
\midrule
\textsc{Blocksworld}\xspace & 119 & 354 & 0.00 & 0.00\\
\textsc{Depots}\xspace      & 200 & 450 & 0.15 & 0.23 \\
\textsc{Driverlog}\xspace   & 183 & 473 & 0.00 & 0.00 \\
\textsc{Logistics}\xspace   & 85 & 446 & 0.12 & 0.01 \\
\ensuremath{\textsc{Satellite}}\xspace   & 117 & 496 & 0.04 & 0.00  \\
\textsc{Zenotravel}\xspace  & 83 & 350 & 0.00 & 0.00 \\
\bottomrule
\end{tabular}
\caption{Hyperparameters of the networks used in our experiments. ${\abs{E}}$ is the size of the vector in output from the embedding layer, ${\abs{LSTM}}$ is the number of neurons in the LSTM layer.}
\label{tab:hyperparameters_sup}
\end{table}

\subsection{Test Set Instances}

Table \ref{tab:test_support_a}  reports the number of GR instances for the test sets used in our experiments.
Domains \textsc{Blocksworld}\xspace and \ensuremath{\textsc{Satellite}}\xspace have fewer instances
 to avoid the overlapping between train and test instances, because of the more limited  space of GR instances in these domains.

\begin{table}[h!]
\centering
\begin{tabular}{lccc}
\toprule
Domain & $\abs{\ensuremath{\textsc{TS}_{\textit{Per}}}\xspace}$ & $\abs{\ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace} $ & $\abs{\ensuremath{\textsc{TS}_{\textit{Rec}}}\xspace}$ \\
\midrule
\textsc{Blocksworld}\xspace & 246 & 3136 & 900 \\
\textsc{Depots}\xspace      & 84 & 6720 & 900 \\
\textsc{Driverlog}\xspace   & 84 & 6720 & 900 \\
\textsc{Logistics}\xspace   & 153 & 6720 & 900 \\
\ensuremath{\textsc{Satellite}}\xspace   & 84 & 5760 & 900 \\
\textsc{Zenotravel}\xspace  & 84 & 6720 & 900 \\
\bottomrule
\end{tabular}
\caption{Number of GR instances per domain for the test sets used in the experiments.}
\label{tab:test_support_a}
\end{table}

\subsection{Additional Material about the Experimental Analysis}\label{detailed_results}

Table \ref{tab:results_extra_our} reports results about the performance of \ensuremath{{\sf GRNet}}\xspace using \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace in terms of Accuracy, Precision, Recall and F-Score.

\vspace{0.2cm}
\noindent
Table \ref{tab:results_extra_per} reports results about the performance of {\ensuremath{\sf LGR}}\xspace using \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace in terms of Accuracy, Precision, Recall and F-Score.

\vspace{0.2cm}
\noindent
Table \ref{tab:results_lama_lpg_appendix} reports results about the performance of \ensuremath{{\sf GRNet}}\xspace in terms of accuracy with the instances in \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace created from plans generated using either \textsc{lpg}\xspace, \textsc{lama}\xspace, or both planners.

\vspace{0.2cm}
\noindent
Table \ref{tab:times} shows the average execution times (milliseconds) of{\ensuremath{\sf LGR}}\xspace and \ensuremath{{\sf GRNet}}\xspace for tests set \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace.

\begin{table}[h!]
\small
    \centering
    \begin{tabular}{ lcccccc }
        \toprule
        Domain & Observations & Plan Length & Accuracy  & Precision & Recall & F1-Score \\
        \midrule
        \multirow{3}{*}{\textsc{Blocksworld}\xspace} & 30\% & 10.1 & 53.4 & 54.4 & 53.2 & 52.1\\
        & 50\% & 17.5 & 69.4 & 71.2 & 70.2 & 69.0\\
        & 70\% & 24.1 & 82.6 & 83.2 & 81.8 & 81.1\\
        \hline
        \multirow{3}{*}{\textsc{Depots}\xspace} & 30\% & 8.4 & 60.9 & 60.9 & 60.8 & 60.6 \\
        & 50\% & 14.6 & 75.8 & 76.6 & 76.4 & 76.3\\
        & 70\% & 20.3 & 86.3 & 87.3 & 87.1 & 87.1\\
        \hline
        \multirow{3}{*}{\textsc{Driverlog}\xspace} & 30\% & 7.0 & 61.4 & 61.7 & 61.9 & 61.4 \\
        & 50\% & 12.2 & 74.8 & 75.9 & 74.1 & 75.0\\
        & 70\% & 17.0 & 83.9 & 84.7 & 84.2 & 83.1\\
        \hline
        \multirow{3}{*}{\textsc{Logistics}\xspace} & 30\% & 7.4 & 65.7 & 67.9 & 66.8 & 66.0 \\
        & 50\% & 12.9 & 78.5 & 79.4 & 78.5 & 78.0 \\
        & 70\% & 17.9 & 87.4 & 87.3 & 86.7 & 86.3\\
        \hline
        \multirow{3}{*}{\ensuremath{\textsc{Satellite}}\xspace} & 30\% & 4.7 & 71.9 & 72.2 & 72.3 & 70.9 \\
        & 50\% & 8.3 & 84.0 & 83.2 & 83.3 & 82.2 \\
        & 70\% & 11.6 & 93.7 & 92.0 & 93.0 & 92.2 \\
        \hline
        \multirow{3}{*}{\textsc{Zenotravel}\xspace} & 30\% & 6.1  & 77.0 & 77.5 & 76.6 & 76.5\\
        & 50\% & 10.6 & 89.7 & 90.0 & 89.5 & 89.5\\
        & 70\% & 14.8 & 96.0 & 96.0 & 95.9 & 95.8\\
        \bottomrule
    \end{tabular}
    \caption{Experimental results about the performance of \ensuremath{{\sf GRNet}}\xspace.
    Columns Observations and Plan Length correspond  to the observability percentage over the full plan and the average number of observed actions in test instances, respectively.
    Each value in column Precision is an average of the precision for the gropus of instances that have the same goal set $\cal G$. Similarly for the values of columns Recall and F1 score.
    }

    \label{tab:results_extra_our}
\end{table}

\newpage
\begin{table}[h!]
    \small
    \centering
    \begin{tabular}{ lcccccc }
        \toprule
        Domain & Observations & Plan Length & Accuracy  & Precision & Recall & F1-Score \\
        \midrule
        \multirow{3}{*}{\textsc{Blocksworld}\xspace} & 30\% & 10.1 & 38.9 & 39.3 & 37.6 & 37.3\\
        & 50\% & 17.5 & 52.9 & 59.6 & 52.4 & 53.0\\
        & 70\% & 24.1 & 76.0 & 79.3 & 75.6 & 76.0\\
        \hline
        \multirow{3}{*}{\textsc{Depots}\xspace} & 30\% & 8.4 & 45.9 & 46.4 & 45.9 & 45.7 \\
        & 50\% & 14.6 & 65.4 & 65.9 & 65.0 & 64.9\\
        & 70\% & 20.3 & 82.1 & 82.3 & 82.1 & 82.1\\
        \hline
        \multirow{3}{*}{\textsc{Driverlog}\xspace} & 30\% & 7.0 & 43.9 & 45.3 & 43.9 & 44.0\\
        & 50\% & 12.2 & 61.5 & 62.7 & 61.6 & 61.7\\
        & 70\% & 17.0 & 79.7 & 80.1 & 79.6 & 79.7\\
        \hline
        \multirow{3}{*}{\textsc{Logistics}\xspace} & 30\% & 7.4 & 50.1 & 51.1 & 50.6 & 50.6 \\
        & 50\% & 12.9 & 67.8 & 68.7 & 68.3 & 68.1 \\
        & 70\% & 17.9 & 82.6 & 83.7 & 83.7 & 83.4\\
        \hline
        \multirow{3}{*}{\ensuremath{\textsc{Satellite}}\xspace} & 30\% & 4.7 & 69.0 & 71.5 & 68.8 & 68.0 \\
        & 50\% & 8.3 & 83.2 & 83.0 & 82.3 & 81.5 \\
        & 70\% & 11.6 & 90.3 & 89.5 & 90.0 & 89.2 \\
        \hline
        \multirow{3}{*}{\textsc{Zenotravel}\xspace} & 30\% & 6.1  & 47.9 & 48.0 & 47.6 & 47.7\\
        & 50\% & 10.6 & 68.8 & 68.9 & 68.9 & 68.9\\
        & 70\% & 14.8 & 87.1 & 87.1 & 87.1 & 87.1\\
        \bottomrule
    \end{tabular}
    \caption{
    Experimental results about the performance of {\ensuremath{\sf LGR}}\xspace.
    Columns Observations and Plan Length correspond  to the observability percentage over the full plan and the average number of observed actions in test instances, respectively.
    Each value in column Precision is an average of the precision for the gropus of instances that have the same goal set $\cal G$. Similarly for the values of columns Recall and F1 score.
    }
    \label{tab:results_extra_per}
\end{table}

\begin{table}[]
    \centering
    \begin{tabular}{ l|ccc|ccc|ccc }
        \toprule
        \multirow{2}{*}{Domain} & \multicolumn{3}{c |}{30\% of  the plan} & \multicolumn{3}{c|}{50\% of the plan} & \multicolumn{3}{c}{70\% of the plan} \\

        & \textsc{lpg}\xspace & \textsc{lama}\xspace & \textsc{lpg}\xspace/\textsc{lama}\xspace & \textsc{lpg}\xspace & \textsc{lama}\xspace & \textsc{lpg}\xspace/\textsc{lama}\xspace & \textsc{lpg}\xspace & \textsc{lama}\xspace & \textsc{lpg}\xspace/\textsc{lama}\xspace \\
        \midrule
        \textsc{Blocksworld}\xspace & 52.4 & 54.4 & 53.4 & 70.0 & 68.8 & 69.4 & 82.3 & 82.9 & 82.6 \\
        \textsc{Depots}\xspace & 61.0 & 60.9 & 60.9 & 76.5 & 75.2 & 75.8 & 87.2 & 85.5 & 86.3 \\
        \textsc{Driverlog}\xspace & 64.9 & 57.9 & 61.4 & 78.1 & 71.5 & 74.8 & 86.2 & 81.6 & 83.9 \\
        \textsc{Logistics}\xspace & 66.8 & 64.6 & 65.7 & 78.5 & 78.4 & 78.5 & 86.8 & 88.0 & 87.4 \\
        \ensuremath{\textsc{Satellite}}\xspace & 72.7 & 71.2 & 71.9 & 84.0 & 83.9 & 84.0 & 93.6 & 93.9 & 93.7 \\
        \textsc{Zenotravel}\xspace & 77.6 & 76.4 & 77.0 & 89.4 & 89.9 & 89.7 & 95.9 & 96.1 & 96.0 \\
        \bottomrule
    \end{tabular}
    \caption{Goal classification accuracy of \ensuremath{{\sf GRNet}}\xspace
    for  \ensuremath{\textsc{TS}_{\textit{PerGen}}}\xspace instances created from plans generated by \textsc{lpg}\xspace, \textsc{lama}\xspace, and both planners over six domains.
    }
    \label{tab:results_lama_lpg_appendix}
\end{table}

\begin{table}[h!]
    \centering
    \begin{tabular}{ l|cc|cc|cc|cc}
        \toprule
        \multirow{2}{*}{Domain} & \multicolumn{2}{c |}{30\%} & \multicolumn{2}{c|}{50\%} & \multicolumn{2}{c |}{70\%} & \multicolumn{2}{c}{Overall mean} \\

        & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace & {\ensuremath{\sf LGR}}\xspace & \ensuremath{{\sf GRNet}}\xspace \\
        \midrule
        \textsc{Blocksworld}\xspace & 832 & 51 & 830 & 50 & 839 & 51 & 834 & 51  \\
        \textsc{Depots}\xspace & 1138 & 65 & 1133 & 82 & 1150 & 67 & 1140 & 71 \\
        \textsc{Driverlog}\xspace & 1171 & 71 & 1170 & 67 & 1163 & 68 & 1168 & 69 \\
        \textsc{Logistics}\xspace & 1263 & 62 & 1261 & 75 & 1262 & 64 & 1262 & 67\\
        \ensuremath{\textsc{Satellite}}\xspace & 1507 & 64 & 1551 & 62 & 1543 & 62 & 1534 & 62\\
        \textsc{Zenotravel}\xspace & 1556 & 55 & 1576 & 52 & 1568 & 53 & 1567 & 54 \\
        \bottomrule
    \end{tabular}
    \caption{Comparison of the execution times of {\ensuremath{\sf LGR}}\xspace and \ensuremath{{\sf GRNet}}\xspace in milliseconds.
    }
    \label{tab:times}
\end{table}