\section{Introduction}
\label{sec-intro}
Automatic License Plate Recognition (ALPR) systems have received much attention to modern transportation services mainly for  automatic management of traffic, parking, toll\hyp station, and other road operations including surveillance and recognition of potential threats~\cite{du2012automatic}.
An ALPR system consists of three main phases: \textit{frame\hyp selection}, \textit{character\hyp segmentation}, and\textit{ optical character recognition} (OCR).
The first phase verifies the existence of any Vehicle License Plate (VLP) in the input frame by extracting any possible character's features from the frame. The second phase separates the characters from the background followed by their recognition in the last phase.
Such a system has more potential for the developing countries or states
requiring higher road\hyp safety measures and better traffic management.\looseness -1

One such potential group of the developing regions is the Bengali speaking  countries and states that still requires a promising ALPR system for the Bengali VLP application.
Unlike English characters, Bengali has more complex features leaving its accurate recognition from a VLP more challenging. Despite being the six largest population  around the world~\cite{ ethnologue2021}, no significant progress can be tracked in the Bengali language countries or states for the ALPR system addressing their inadequate road\hyp safety measures and poor traffic management. Besides, the performance of the prominent ALPR models developed for the English VLP is also unknown for the Bengali VLP recognition application.\looseness -1

Recent ALPR systems for English VLP captures the employment of promising machine learning models. Such models mainly detects the VLP in the given image (or a video frame) followed by the recognition of the text information on the plate.
A variety of models are stemmed from the need for improving the classification accuracy, robustness to environmental artifacts, and computational efficiency.
For example, the Convolutional Neural Network (CNN) based bounding box detectors were developed with regression algorithms~\cite{bulan2017segmentation, Xueccv2018, laroca2018, wang2019, silva2021}, manual annotation~\cite{Zhuang_2018_ECCV}, and transfer learning~\cite{c12}.
Other development with  recurrent neural network (RNN) based architectures include the  Bidirectional Long short-term memory (BiLSTM) based models~\cite{zou2020robust}, and real\hyp time object detection algorithm, YOLO (You Only Look Once) and its variants based models~\cite{chen2019automatic, al2019efficient, yolo}. The end-to-end cascaded and unified architectures of the neural networks were also studied for the ALPR systems~\cite{Hsu2013, Montazzolli2017, Li2019, onim2020traffic}.

However, Unlike the English VLP, the effort in developing the ALPR systems for the Bengali VLP is particularly limited.
Despite much interest in developing Bengali handwriting, scripts and character recognition in general, it has not captured the ALPR system yet.
A few notable developments of the ALPR systems for the Bengali VLP include the feature extraction based on the digital curvelet transform~\cite{majumdar2007bangla},
Tesseract OCR~\cite{hasnat2009integrating}, and CNN with  Adam optimizer~\cite{rabby2018bornonet}.

In summary, no Deep Neural Networks (DNN) model can be tracked in the literature that can detect the Bengali VLP and recognize its characters simultaneously.
A few models only focused on the Bengali characters recognition~\cite{ hasnat2009integrating, rabby2018bornonet}, but their performance is unknown for ALPR system.
Although Tesseract and BornoNet have relatively high character recognition accuracy, they are not suitable for real\hyp time applications like ALPR due to higher processing time.
Unlike the above models, Onim~\emph{et al}\onedot\cite{onim2020traffic} recently combined YoloV4 for detection of VLP, and Tesseract as OCR engine.
However, the model requires reasonably higher time for number plate detection and character recognition.
All these mean that the models  developed for either Bengali VLP detection or Bengali characters recognition are generally limited with \textit{low detection and recognition accuracy}, and \textit{high computational complexity} making them unsuitable for a real~\hyp time application.

In this paper, we, therefore, report an ALPR system with a new DNN model that we call Bengali License Plate Network (BLPnet) (Sec.~\ref{sec-our ALPR}).
BLPnet is constructed to have three primary phases.
Particularly, the contributions with the proposed three\hyp phase ALPR system can be summarized as follows:

\begin{itemize}
    \item Building on the NASNet\hyp Mobile backbone architecture, the first phase employs more \textit{dense} and \textit{pulling} layers on the network\hyp head to identify the vehicles more efficiently with a region of interest bounding\hyp box (Sec.~\ref{subsec-boundingbox}).

    \item The second phase includes an InceptionV3 architecture customized with several \textit{dense} and \textit{pulling} layers to detect the VLP in the bounding\hyp box region (Sec.~\ref{subsec-VLPdetect}).

    \item The third phase introduces a new Bengali\hyp OCR engine (Sec.~\ref{subsec-OCR}).

    Unlike the existing ALPR systems that cannot adequately tackle the artifacts like motion\hyp blur and non\hyp uniform shadow, the proposed engine  employs de\hyp blurring and Region Scalable Fitting (RSF) based segmentation, which are optionally invoked (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot, when characters are not recognized) to better tackle the intensity inhomogeneity.
    \item Additionally, detecting vehicle regions (first phase) prior to VLP (second phase), as considered in the proposed ALPR system, would significantly reduce computational cost and false\hyp positives making the system faster and more accurate.
\end{itemize}

\section{Related ALPR Systems}
Development of several prominent ALPR systems for English VLP can be tracked in the literature so far.
Getting the VLP information requires the localization and detection of the VLP. In CNN based approach, features are extracted for VLP and then object localization is done with bounding box detectors and regression algorithms~\cite{Xueccv2018, wang2019}.
Silva~\emph{et al}\onedot~\cite{silva2021} and Laroca~\emph{et al}\onedot~\cite{laroca2018} used bounding box regression to localize the object of interest with image coordinates using YOLO. These are One-stage detection networks and relatively faster than other detectors. But these approach is computationally inefficient as it requires training \textit{darknet} backbone of over $27M$ parameters.

Alternatively, VLP can also be semantically segmented using a either CNN or deep segmentation network.
Bulan~\emph{et al}\onedot\cite{bulan2017segmentation} proposed a real\hyp time complete CNN model  having higher adaptability to tackle the scene\hyp variations  requiring minimum possible human\hyp intervention.
An image\hyp based classification was proposed with weak\hyp Snow and strong\hyp CNN classifiers.
The first classifier was initially used for the region\hyp tagging of the possible number\hyp plates, and the strong\hyp CNN refined the classification outputs for the actual number\hyp plates
Zhuang~\emph{et al}\onedot~\cite{Zhuang_2018_ECCV} proposed such method where they segmented the VLP and the characters in it. This approach also needs a manual annotation of a characters and LP. Besides these segmentation based networks, segmentation free Networks have also been proposed.
Wang~\emph{et al}\onedot~\cite{wang2019} proposed a convolutional Recurrent Neural Network (CRNN) and proposed a multi-task license plate detection and recognition (MTLPR) model.
Zou~\emph{et al}\onedot~\cite{zou2020robust} used Bi-LSTM and directly localized characters without segmentation.\looseness -1

Both end-to-end cascaded and unified architectures have been used for ALPR system. Li~\emph{et al}\onedot~\cite{Li2019}
used end-to-end trainable LP detection model that remains error-prone to the challenging conditions and also depends heavily on architectural design.
On the other hand, cascaded architectures have the flexibility of tuning and testing. Hsu~\emph{et al}\onedot~\cite{Hsu2013} and Montazzolli~\emph{et al}\onedot~\cite{Montazzolli2017} proposed cascaded architecture for real time VLP detection and recognition.

Huang~\emph{et al}\onedot~\cite{c12} introduced Transfer learning into their CNN model to tackle the limited labeled data for a high detection accuracy of VLP of  93\%.
Later, with a modified YOLOv3 model with 10 CNNs, Chen~\emph{et al}\onedot\cite{chen2019automatic} attempted to detect and recognize the VLP  with an accuracy of about 98\% and  78\%, respectively in a variety of conditions (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot, rainy backgrounds, darkness and dimness, and varied colors and saturation of photos).
In addressing several other environmental artifacts like
poor\hyp contrast and noisy images resulting from foggy, distortion and dusty conditions,
Al-Shemarry~\emph{et al}\onedot\cite{al2019efficient} improved the accuracy of VLP detection using a  binary descriptor with contrast enhancement, leaving their proposed system limited to real\hyp time application with poor frame\hyp rate and higher computational cost.
To tackle the motion artifacts (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot, blurry input images),
Zou~\emph{et al}\onedot\cite{zou2020robust} recently proposed a Bi-LSTM model with  robust blur\hyp kernel estimation that offered a reasonably lower accuracy of 79.55\%.
Another recent development with object bounding\hyp box detection for ALPR system is YOLOv4 model~\cite{yolo}, which is
a deep convolutional network that can localize VLPs with immense frame rate and high accuracy.\looseness -1

Besides, on Bengali character detection and recognition,
Majumdar~\emph{et al}\onedot~\cite{majumdar2007bangla} proposed a feature extraction based on the digital curvelet transform.
Separate $k$\hyp nearest neighbor classifiers were trained using the curvelet coefficients of an original image and its morphologically changed copies.
The overall recognition accuracy of 96.8\% was reported
while trained for twenty popular Bengali fonts and tested for different font sizes.
Later, Hasnat~\emph{et al}\onedot~\cite{hasnat2009integrating} attempted to successfully introduce Bengali script in Tesseract OCR.
Tesseract can learn new languages or alphabets with training along with existing languages for character recognition.
Recently, Rabby~\emph{et al}\onedot~\cite{rabby2018bornonet} presented a new 13-layer CNN model called BornoNet for Bengali characters recognition. The model was designed with two sub-layers optimized using adam optimizer
Despite having relatively higher accuracy in the models for Bengali characters recognition, their applications to an ALPR system is not yet studied.
In addressing the potential gap of having a promising ALPR system for Bengali VLP, we therefore develop and present real\hyp time End\hyp to\hyp End  DNN based model that we call BLPnet.

\section{A New ALPR System with BLPnet}
\label{sec-our ALPR}
This section presents a newly developed BLPnet model for ALPR system.
As illustrated in Fig.  \ref{Overall block diaram}, the model consists of three main phases that are cascaded together to detect vehicle, VLP and recognize the characters in it.
It takes real-time video frames as input consisting of vehicles and its surroundings.
In the first phase, the vehicle region is detected and distinguished using a bounding\hyp box
as the smallest possible rectangle with vertical and horizontal sides that entirely surround the vehicle captured in a frame.
The frames containing the bounding boxes are fed to the VLP detection model along with the coordinates of the bounding\hyp box.
Later, the second phase detects and extracts the VLP contour
followed by the character  localization and recognition in the final phase.
Consideration of intensity\hyp inhomogeneity resilient segmentation and effective deployment of these cascaded modules would reduce the computational cost and false\hyp positive predictions improving the accuracy and reliability of the system for a real\hyp time application.
These phases are now discussed below with more technical details.

\begin{figure*}[htp]
\centering
\includegraphics[width=0.95\columnwidth]{Fig-1}
\caption{Key processes of the proposed ALPR system}
\label{Overall block diaram}
\end{figure*}

\subsection{An Extended Vehicle Bounding\hyp Box Detector}
\label{subsec-boundingbox}
A NASNet\hyp Mobile is extended and used as the backbone of our DNN followed by six hidden layers to improve the vehicle detection accuracy (see Fig.~\ref{fig-BBox}).
The NASNet\hyp Mobile is more computationally efficient than its counterpart like the ResNet-based architecture that needs higher hardware configuration, making it unsuitable for large\hyp scale deployment~\cite{c13}.
Thus, our model can be implemented on Field-Programmable Gate Array (FPGA) and other embedded devices to balance between hardware resources and processing speed. The TensorFlow (TF) v2.4 framework was used for the overall construction of our CNN model.
The backbone was pre-trained with the Imagenet dataset \cite{imagenet}.
Three hidden layers added as network head were fully connected layers followed by a dropout layer. A few important processing and considerations are briefly discussed below.\looseness -1

\begin{figure*}[htb]
\centering
\includegraphics[width=0.6\columnwidth]{Fig-2}
\caption{Vehicle Bounding\hyp box detector of the proposed ALPR system}
\label{fig-BBox}
\end{figure*}

\paragraph{Dataset Collection}
For the training of vehicle bounding\hyp box detection model, we used the Cars dataset~\cite{stanford} of Stanford AI Lab.
There are 16,185 images in the dataset, representing 196 different car classifications. The data is divided into 8,144 training images and 8,041 testing images, with about a 50-50 ratio between the two classes. Typically, classes are organized by Make, Model, and Year, such as 2012\hyp Tesla\hyp Model\hyp S or 2012\hyp BMW\hyp M3\hyp coupe. Each image in the dataset has also  a  bounding\hyp box output value as an ideal reference that suits our requirement.

\paragraph{Pre-Processing}
We considered a total of  8,041 training images and passed them through two preprocessing stages. At first image augmentation was done to simulate challenging conditions which the real\hyp time traffic video might face. The applied augmentations include: rotate, shift, flip, varying contrast, blur, and salt and pepper noise.
After augmentation, label-encoding and  training-validation split are done. In this case, we made an 80-20 split. After the splitting, the model is trained and tuned.

\paragraph{Transfer Learning (TL)}
A pre\hyp trained model (NASNet-Mobile) as the backbone architecture of the proposed model ensures optimised training.
Unlike the isolated learning paradigm, Transfer learning allows extending the pre\hyp learned features to adequately address the new related challenges.
Particularly, an average pooling layer usually can extract the pretrained models edge detection capability.
From the final layer of shape $(10\times10\times1056)$, an average pooling layer is thus added followed by a group of fully connected dense layers and dropout layers.
 Finally the output layers identifies the object and generates the bounding box coordinates.
Thus, our model was pre\hyp trained with over fourteen  million images from the ImageNet dataset, which is now capable of categorising images into over 1000 different classes with transfer learning ability.
As a result, it has accumulated a library of rich feature sets for a wide variety of classes.

\paragraph{Hyper-parameters and Tuning}
The training hyper-parameters are given in Table~\ref{tab:training_params}.
Finally, the parameters with the best accuracy have been chosen after several trials and errors.
The optimizer that suits best for the VLP detection is Stochastic gradient descent(SGD).
We have also tweaked several parameters of the optimiser
for our dataset and defined early stopping to reduce the training time.
Additionally, Reduced Learning Rate functionality has been used to further expedite the training process.

\begin{table}[htb]
    \centering
    \caption{Comparison of the training hyper\hyp parameters}
    \setlength{\tabcolsep}{4pt}
    \resizebox{\linewidth}{!}{
    \begin{tabular}{l|cc|cc|c|c|c}
    \toprule
       \bf Model	& \multicolumn{2}{c|}{\bf Ours (BLPnet)}	& \multicolumn{2}{c|}{\bf Onim~\emph{et al}\onedot~\cite{onim2020traffic}}	& \bf Zou~\emph{et al}\onedot~\cite{zou2020robust}	& \bf \makecell{Hendry \&\\ Chen~\cite{chen2019automatic}}	& \bf Li~\emph{et al}\onedot~\cite{LI201814}\\
        \midrule
        \bf Parameters &	Vehicle	& VLP	& Vehicle &	VLP	& VLP	& VLP	& VLP\\
        \midrule
        \bf Backbone	& \makecell[l]{NASNet\\Mobile}	& InceptionV3	& Darknet	& Darknet	& MobileNetV3	& Darknet & - \\
        \midrule
        \bf Classes &	1 &	1 &	4 &	1 &	37 &	60 &	37\\
        \midrule
        \bf Batch size	& 32 &	64 &	32 &	64 &	- &	- &	32\\
        \midrule
        \bf Image Shape &	$300\times300$ &	$200\times200$ &	$224\times224$ &	$128\times128$ &	$152\times56$ &	$256\times256$ &	$24\times24$ \\
        \midrule
        \bf Optimizer &	SGD &	Adam &	Adam &	Adam &
        SGD &	Adam &	Adam\\
        \midrule
        \bf Epochs &	300 &	300 &	100 &	3500 & - &		600 & - \\
    \bottomrule
    \end{tabular}
    \label{tab:training_params}
    }
\end{table}

\subsection{VLP Detection with Inception-v3}
\label{subsec-VLPdetect}
We have also extended the
Inception-v3 DNN architecture and employed it as the backbone of the second phase process of BLPnet for VLP detection as illustrated in Fig.~\ref{lp-BBox}.
Inception-v3 is a CNN designed from the Inception family that transports label\hyp  information to a lower level of the network using label\hyp smoothing, factorised $7\times 7$ convolutions, an auxiliary classifier, and batch\hyp  normalization for layers in the side head.
Inception architecture being 42 layers deep has 2.5 times lower  computational cost than that of GoogLeNet and can function well even when memory and computational budgets were limited~\cite{szegedy2016rethinking}.

We have used an annotated dataset of 1500 training and 300 validation images to train our model.
From real\hyp world video footage taken on roadways in Dhaka, Bangladesh, with varying road conditions, including high traffic congestion, the VLP were detected with reasonably higher accuracy (see Sec.~\ref{sec-result}).
The training began with a number of iterations set at 6000.
The average loss did not diminish significantly after 2000 iterations,
and after each 1000 cycles, the training was capable of taking weight backup.
To reduce training duration, we thus used early stopping.

\begin{figure*}[htb]
\centering
\includegraphics[width=0.6\columnwidth]{Fig-3}
\caption{License Plate Bounding\hyp box detector of the proposed ALPR system}
\label{lp-BBox}
\end{figure*}

\subsection{A New Bengali\hyp OCR Engine}
\label{subsec-OCR}
This step recognised characters in the extracted VLP from the previous step.
Similar to~\cite{9003211}, we addressed the character recognition as an object recognition problem.
By treating characters as objects, the character segmentation and recognition steps are integrated together that consists of a total 5 convolutional 2D layers with 16, 32, 64, 128 and 512 nodes.
Each of them includes a preceding max-pooling layer with a pool size of 2 and a dropout layer with a drop rate of 20\%.
The dropout layers were used to prevent the model from over-fitting.
The kernel size of 2 is considered with necessary padding.
Relu activation function is used for all the backbone layers.
Finally, a global\hyp average pooling layer is included before the output layer followed by two dense layers.
Two important considerations of this model are the conditional use of de\hyp blurring filter and intensity inhomogeneity invariant segmentation of the character.

Blurring is known to be a key challenge in VLP recognition resulting from wrong focusing or vehicle's motion that requires
restoration of the blurred images
for a higher accuracy in character recognition.
A set of filters like  Wienner filter, Total Variation (TV) deblurring, Haar deblurring and a combination of both is considered  for restoring any blurred image to its estimated original version.
For an input image,  $Im[n]$, let the output of the filter is  $x[n]$ as in~\eqref{eqn_xn}, where $G(z)$ is the filter function, $N$ is the number of past taps with coefficient $a$, $s[n]$  is the reference signal, and $e[n]$ is the residual error.
The criteria to minimize the mean squared error (MSE) for the coefficient is defined in~\eqref{eq-a}, where $E[.]$ is the expectation operator.
On the other hand, Haar deblurring~\cite{DBLP:journals/corr/abs-1905-01003} uses a deconvolution filter and reduce the MSE as in~\eqref{eq-har}, where  $I(i,j)$ is the expected and $G(i,j)$ is the desired deblurred images.
In our model, we have designed a conditional deblurring process invoked when not enough number of characters are detected in the OCR.
An iterative approach is employed for this process with a varying threshold
based on the Fast-Iterative Shrinkage Thresholding Algorithm (FISTA)~\cite{doi:10.1137/080716542} in the VLP image.\looseness -1

\begin{subequations}
\begin{equation}
\label{eqn_xn}
x[n]=\sum_{i=0}^{N} a_{i} \times \operatorname{Im}[n-1]
\end{equation}
\begin{equation}\label{eq-a}
a_{i}=\operatorname{argmin}\left(E\left[e^{2}[n]\right]\right)\\
\end{equation}
\end{subequations}\vspace{2pt}

\begin{equation}\label{eq-har}
M S E=\frac{1}{m \times n} \sum_{i=0}^{m-1} \sum_{j=0}^{n-1}[I(I, J)-G(i, j)]^{2}
\end{equation}

Additionally, we have considered the intensity inhomogeinity invariant segmentation of the characters. Thus, we extracted the characters from the VLP with two segmentation models: Chan~\&~Vese (CV) model and Region-Scalable Fitting (RSF) model (see Table~\ref{Table for Performance Analysis of OCR} in Sec.~\ref{sec-result}). As most of the images for VLP have very high contrast, CV model performed well in approximating the intensities of the foreground or object and background.
In case of intensity inhomogeneity, CV model failed to segment the characters, where RSF model was used.

The complete architecture is illustrated in Fig.~\ref{OCR}.
The output layer contains 60 classes to capture the total number of Bengali alphabets and digits.
The output layer uses \textit{Softmax} as an activation function.
We have also used CMATERdb version 3.1.2 for training the alphabets for the digits. A total alphabet\hyp images of  15000 and digits of 500 have been split as follows: 9900 images for training,   2600 images for validation, and 3000 images for testing.

\begin{figure*}
    \centering
    \includegraphics[width=0.8\linewidth]{Fig-4}
    \caption{The newly developed OCR engine for BLPnet}\vspace{-3pt}
    \label{OCR}
\end{figure*}

The images are then augmented so that the model can identify images on non-ideal conditions. The following augmentations are considered: (\textit{i})~10\% shift in left and 10\% shift in right; (\textit{ii})~10\% shift in vertical and 10\% shift in horizontal; (\textit{iii})~7.5$\degree$ rotation in left and 7.5$\degree$ in right; (\textit{iv})~zoom range 20\%; and, (\textit{v})~Shear range 20\%.
Besides, each character illustrated in Fig.~\ref{Character Class} is resized to $16\times 16$ after segmentation before passing through the network.
As a result,  this network can work with variable size characters.

\begin{figure}[thb]
\centerline{\includegraphics[width=0.55\linewidth]{Fig-5}}
\caption{Bengali character classes for the OCR training}\vspace{-3pt}
\label{Character Class}
\end{figure}

\section{Results and Analysis}
\label{sec-result}
We now present the performance of the proposed BLPnet model for the ALPR system, and compare its performance with some related models:
Li~\emph{et al}\onedot~\cite{LI201814}, Hendry \& Chen~\cite{chen2019automatic}, Zou~\emph{et al}\onedot~\cite{zou2020robust}, and Onim~\emph{et al}\onedot~\cite{onim2020traffic}.
This section starts with the  analysis of total learnable parameters of the model followed by the
 results discussed in three main parts, namely, for vehicle detection, VLP detection and finally, for the OCR and word\hyp mapping.

\subsection{Total Learnable Parameters}
A dense layer with \textit{m} input nodes and \textit{n} output nodes will have a total of $(n+1)\times m$ learnable parameters.
The pooling layers and dropout layers does not learn anything. For convolutional layers with \textit{p} feature maps in input and \textit{q} feature maps as output having a filter size of $i\times j$ will have a total learnable parameters of $(i\times j\times p + 1)\times q$.
 The distribution of trainable parameters for bounding box detector across the network is noted in Table.~\ref{tab:bbox_param}.
 The total number of trainable parameters is calculated as shown in Equation.~\eqref{param}.
 Similarly, Table.~\ref{tab:ocr_params} shows the trainable parameters for our OCR engine. Here \textit{Layers} denotes the name and operation of added layers, \textit{Shape} denotes the input shape of tensors for that particular layer and finally the trainable parameter for that layer.

 \begin{table}[!th]
     \centering
     \caption{Trainable parameters for bounding box detector}
     \renewcommand{\arraystretch}{1.05}
	\setlength{\tabcolsep}{10pt}
	\resizebox{0.7\linewidth}{!}{
     \begin{tabular}{llc}
\toprule
\bf Layers &\bf Shape &\bf \makecell[c]{Trainable\\ Parameters} \\
\midrule
Global\_AveragePooling2D & (None, 1056) & 0 \\
Dense & (None, 256) & 270592 \\
Dense & (None, 256) & 270592 \\
Dense & (None, 128) & 32896 \\
Dropout & (None, 256) & 0 \\
Dropout & (None, 128) & 0 \\
Dense & (None, 128) & 32896 \\
Dense & (None, 64) & 8256 \\
Dense & (None, 64) & 8256 \\
Dropout & (None, 64) & 0 \\
Dropout & (None, 64) & 0 \\
Dense & (None, 32) & 2080 \\
Dense & (None, 32) & 2080 \\
Dense & (None, 2) & 66 \\
Dense & (None, 4) & 132 \\
\midrule
\bf Total & & 627k\\
\bottomrule
\end{tabular}
}
     \label{tab:bbox_param}
 \end{table}

\begin{equation}
\begin{aligned}
    N_{detector} &=2\times((1056+1) \times 256+(256+1) \times 128+(128+1) \times 64 \\
    & +(64+1) \times 32+(32+1) \times 2) \\
    &=627k
\end{aligned}
\label{param}
\end{equation}

\begin{table}[!t]
    \centering
    \caption{Trainable parameters for OCR engine}
\renewcommand{\arraystretch}{1.05}
	\setlength{\tabcolsep}{8pt}
	\resizebox{0.6\linewidth}{!}{
\begin{tabular}{llc}
\toprule
\bf Layers &  \bf Shape & \bf \makecell[c]{Trainable\\ Parameters}  \\
\midrule
Conv2D  & (None, 63, 63, 16) & 80 \\
MaxPooling2D & (None, 31, 31, 16) & 0 \\
Dropout & (None, 31, 31, 16) & 0 \\
Conv2D  & (None, 30, 30, 32) & 2080 \\
MaxPooling2D & (None, 15, 15, 32) & 0 \\
Dropout & (None, 15, 15, 32) & 0 \\
Conv2D  & (None, 14, 14, 64) & 8256 \\
MaxPooling2D & (None, 7, 7, 64) & 0 \\
Dropout & (None, 7, 7,64)  & 0 \\
Conv2D  & (None, 6, 6, 128) & 32896 \\
MaxPooling 2D & (None, 3, 3, 128)  & 0 \\
Dropout  & (None, 3, 3, 128)  & 0 \\
Conv2D  & (None, 2, 2, 256)  & 131328 \\
MaxPooling 2D  & (None, 1, 1, 256)  & 0 \\
Dropout  & (None, 1, 1, 256)  & 0 \\
Flatten  & (None, 256) & 0 \\
Dense & (None, 256)  & 65792 \\
Dense & (None, 512)  & 131584 \\
Dropout & (None, 512) & 0 \\
Dense  & (None, 60) & 25650 \\
\midrule
\bf Total &  & 397K\\
\bottomrule
\end{tabular}
}
    \label{tab:ocr_params}
\end{table}

\subsection{Vehicle Detection}
We trained the vehicle bounding\hyp box model for vehicle detection for 1000 epochs with early stopping employed to stop the model training if the accuracy improvement is negligible.
The hardware specification to train the model is: Tesla K80 2496 CUDA cores GPU with 12~GB GDDR5 VRAM, Processor @2.3Ghz (4 core, 8 threads), and RAM of 48~GB.

\begin{figure}[!th]
\centerline{\includegraphics[width=.8\linewidth]{Fig-6}}
\caption{Loss vs Epoch curve for the first 50 epochs}
\label{Loss vs Epoch curve}
\end{figure}

\begin{figure*}[htbp]
\centering
    \subfloat[Detected vehicles (marked with green bounding\hyp boxes) ]{%
    \includegraphics[width=0.9\columnwidth, height = 3 in]{Fig-7a}\label{Bounding box output}}
    \\
    \subfloat[VLP detection (marked with blue\hyp colour bounding box)]{%
    \includegraphics[width=0.9\columnwidth, height = 3in]{Fig-7b}\label{cars}}
\caption{Example of detected vehicles and VLP from real\hyp time captured video footage with different orientations of the vehicles}
\end{figure*}

We chose the criteria for our evaluation matrices of bounding\hyp box coordinates to be MSE and accuracy.
From the loss vs epoch curve in Fig.~\ref{Loss vs Epoch curve}, we observe that the model can quickly converge to a lower loss value. After a while, the improvement gets very slow.
The best training loss and MSE of bounding\hyp box were 0.0130 and 0.0155, respectively and the validation loss and MSE were 0.0148 and 0.0152, respectively.
Overall training parameters with their respective epochs is presented in Table~\ref{table of vb training}.
We observe in this table that after 300 epochs, the rate of change in loss and accuracy is becoming stable.
Compared to the YOLO model, the hyper-parameters were kept similar as shown earlier in Table~\ref{tab:training_params}.
All these observations suggest that the our proposed model would effectively detect the vehicles in the input video clips. The final predicted output of the model applied in random vehicles is shown in Fig.~\ref{Bounding box output} with successful demarcation of green boxes.

\begin{table*}[!th]
\caption{Performance of vehicle detection with bounding\hyp box}
\centering
\renewcommand{\arraystretch}{1.05}
	\setlength{\tabcolsep}{3pt}
	\resizebox{\linewidth}{!}{
\begin{tabular}{cccccccccc}
\toprule
\multirow{3}{*}{\textbf{Epochs}} &
\multirow{3}{*}{\textbf{Loss}} & \textbf{Box} & \textbf{Class} & \textbf{Box} & \textbf{Class} & \multirow{3}{*}{\makecell{\bf Validation \\ \bf Loss}} & \makecell{\bf Validation\\ \bf box} & \makecell{\bf Validation\\\bf  class} &  \makecell{\bf Validation\\ \bf box} \\
 &  & \textbf{Output} & \textbf{Output} & \textbf{Output} & \textbf{Output} &  & \textbf{Output} & \textbf{Output} & \textbf{Output} \\
  &  & \textbf{Loss} & \textbf{Loss} & \textbf{MSE} & \textbf{Accuracy} &  & \textbf{Loss} & \textbf{Loss} & \textbf{MSE} \\
    \midrule

50&
0.0056 &
0.0176 &
0.0059 &
0.0176 &
75.59 &
0.0057 &
0.0175 &
0.0057 &
0.0175 \\

100&
0.0056 &
0.0171 &
0.0058 &
0.0171&
85.89 &
0.0056 &
0.0140 &
0.0057 &
0.0140 \\

150&
0.0056 &
0.0169 &
0.0058 &
0.0169 &
89.19 &
0.0055 &
0.0138 &
0.0056 &
0.0138 \\

200&
0.00569 &
0.0173 &
0.0057 &
0.0173 &
90.55 &
0.0055 &
0.0175 &
0.0055 &
0.0175 \\

250&
0.00569 &
0.0176 &
0.0056 &
0.0176 &
94.74 &
0.0054 &
0.0175 &
0.0055 &
0.0175 \\

300&
0.00569 &
0.0176 &
0.0053 &
0.0176 &
96.99 &
0.0055 &
0.0175 &
0.0054 &
0.0175 \\

350&
0.00569 &
0.0176 &
0.0052 &
0.0176 &
97.02 &
0.0055 &
0.0175 &
0.0054 &
0.0175 \\

400&

0.00582 &
0.0130 &
0.0050 &
0.0155 &
97.05 &
0.0056 &
0.0148 &
0.0058 &
0.0152 \\
\bottomrule
\end{tabular}
}
\label{table of vb training}
\end{table*}

\subsection{VLP Detection}

Our model was evaluated using both real\hyp time and pre-recorded video clips. Table \ref{tab:training_params} shows the hyper\hyp parameters used to train the network. Our model did not over\hyp fit due to the use of adequate \textit{dropout} and \textit{pooling} layers, and thus, it took little time to train and converge to an optimum detection level. During the evaluation, the performance of the model was monitored in real\hyp time.

Our algorithm successfully detected VLP as a few examples are illustrated in Fig.~\ref{cars}, where VLP is marked with blue bounding box.
A VLP was detected from varying orientation of the detected vehicle and at that time, 17 frames per second processing speed was maintained on average.
Such accurate detection was observed for all the testing video clips with resolution $1920\times  1080$ and  frames per second (\textit{fps}) is of 15 and 20.\looseness -1

Similar to the vehicle detection with  bounding box (Sec.~\ref{subsec-boundingbox}), the VLP detection training loss and MSE of bounding\hyp box were 0.013 and 0.016, respectively and the validation loss and MSE were 0.015 and 0.015, respectively.
No false\hyp positive detection was recorded during the process which justifies the effectiveness of previous phase.  \looseness-1

\subsection{Character Recognition}
The overall procedure can be broken down into three steps. First the preprocessing step, where the image is filtered based on the sharpness.
Later in character segmentation, the fast active contour model segments the characters as objects. Finally, the character recognition followed by a word mapping. Here, for each  group of key characters recognised, a predefined word is mapped for them. Our model's generated character output is shown in Fig.~\ref{Example results OCR}.

\begin{figure}[!t]
    \centering
 \vspace{-0.3cm}
    \includegraphics[width=0.985\linewidth]{Fig-8}
    \caption{Example of input and outputs of the newly developed OCR engine (From \textit{top}, input images in \textit{first row}, CV\hyp segmented images in \textit{second row},  LSF\hyp segmented images in \textit{third row}, and OCRs and word\hyp mapping in \textit{fourth row})}
    \label{Example results OCR}
\end{figure}

\begin{table}[!t]
\caption{OCR performance for the CV and RSF based segmented characters}
\centering
\renewcommand{\arraystretch}{1}
	\setlength{\tabcolsep}{12pt}
	\resizebox{0.95\columnwidth}{!}{
    \begin{tabular}{ccccc}

\toprule
\textbf{Segmentation} & \textbf{No of} &
\textbf{Accuracy} &
\textbf{Time taken} &
\textbf{Time taken}\\
\textbf{Model} & \textbf{characters} &
\textbf{of OCR} &
\textbf{for OCR} &
\textbf{for Tesseract}\\
& \textbf{extracted} & \textbf{(in \%)} & \textbf{(in seconds)} & \textbf{(in seconds)}
\\
    \midrule
\textbf{CV model} & 4 & 90 & 0.302 & 0.402\\
& 5 & 83 & 0.395 & 0.548\\
& 6 & 81 & 0.432 & 0.701\\
& 8 & 80 & 0.502 & 0.705\\ \midrule
\textbf{RSF model} & 4 & 95 & 0.256 & 0.402\\
& 5 & 93 & 0.312 & 0.548\\
& 6 & 90 & 0.333 & 0.701\\
& 8 & 89 & 0.398 & 0.705\\\bottomrule

    \end{tabular}
    }

    \label{Table for Performance Analysis of OCR}
\end{table}

\begin{table}[!t]
    \caption{Average testing accuracy in characters recognition applications}
    \centering
    \renewcommand{\arraystretch}{1.2}
    \setlength{\tabcolsep}{10pt}
    \resizebox{0.95\linewidth}{!}{%
    \centering
    \begin{tabular}{lll}
        \toprule \multicolumn{1}{l}{ \bf Model } & \bf Accuracy & \bf Recognising object \\
        \midrule

        Hendry \& Chen~\cite{chen2019automatic} & $78.2\%$ &  VLP \& characters (English)\\
        Zou~\emph{et al}\onedot~\cite{zou2020robust} & $79.5\%$ & Blurry text recognition (English)\\
        Onim~\emph{et al}\onedot~\cite{onim2020traffic} & $90.51\%$ & VLP (Bengali)\\

        \textbf{BLPnet (proposed)} & $\mathbf{95.0\%}$ & Real\hyp time VLP \& characters (Bengali)\\
        \bottomrule
    \end{tabular}
    }
    \label{tab_comp}
\end{table}

\begin{table}[!t]
    \centering
    \caption{Comparison of features of the ALPR systems}
    \setlength{\tabcolsep}{8pt}
	\resizebox{\linewidth}{!}{
    \begin{tabular}{lp{3cm}cccp{3cm}}
\toprule
\bf Model  &
\bf  \makecell[c]{VLP\\False--Positive\\Detection} &
\bf \makecell{Character \\rotational \\ invariant?} &
\makecell{\bf Trainable \\\bf parameters\\ \bf ({\it \textbf{Million}})} & \bf \makecell{Processing \\time\\ (\textit{s})}&
\bf Hardware\\
\midrule
\makecell[l]{Hendry \&\\ chen~\cite{chen2019automatic}} & Unknown & no & 8.8 & 0.8 -- 1.0 & \makecell[l]{Core i7\\ Nvidia GTX\\ 970 4GB GPU}\\
\midrule
Li~\emph{et al}\onedot~\cite{LI201814} &
\makecell[l]{Heuristically \\minimized \\with CNN} & no & 1.0 &
2.0 -- 3.0 & \makecell[l]{Core i5\\ Nvidia Tesla\\ k40c 4GB GPU}\\
\midrule
Zou~\emph{et al}\onedot~\cite{zou2020robust} & Unknown & yes & 1.9 & Unknown & \makecell[l]{Core i5\\ Nvidia Titan\\ 12GB GPU}\\
\midrule
Onim~\emph{et al}\onedot~\cite{onim2020traffic} & \makecell[l]{Eliminated using\\ 2\hyp phase detection} & no & 27& 0.7 -- 1.0 &\makecell[l]{Core i5\\ Nvidia Tesla\\ T4 6GB GPU} \\
\midrule
\bf \makecell[l]{BLPnet \\(Ours) }& \bf \makecell[l]{Eliminated using\\ 2\hyp phase detection}  & \bf yes & \bf  0.97 &
\bf 0.32 -- 0.52  & \bf  \makecell[l]{Core i5\\ Nvidia Tesla \\k80 6GB GPU}\\
\bottomrule
\end{tabular}
    \label{tab:major_comp}
    }
\end{table}

Our proposed OCR engine demonstrated reasonably better accuracy  as it treated the characters as objects. The accuracy varies from 80\% to 95\% based on the number of characters extracted. The time it takes to extract a character varies between 0.302 to 0.502 seconds. The model is weak against some vowel pairs. This weakness was tackled successfully with word mapping (see Fig.~\ref{Example results OCR}). Detailed performance analysis is given in Table~\ref{Table for Performance Analysis of OCR}. The accuracy decreases with the increase in number of characters.
A comparison with related OCR engines is shown in Table~\ref{tab_comp}.
The other performance of the proposed model is also compared  with the other prominent models in Table~\ref{tab:major_comp}.
Considering the nullified False\hyp Positives, lower computational complexity, reduced trainable parameters, and reasonably higher accuracy, the proposed model has demonstrated a higher potential for a real\hyp time ALPR application.

 \section{Conclusion}
 The development of an ALPR system is a timely requirement for modern transportation services.
 However, despite being the sixth largest population  around the world~\cite{ ethnologue2021}, no significant progress can be tracked in the Bengali language countries or states for the ALPR system addressing their inadequate road\hyp safety measures and poor traffic management.
 To this end, we have presented a  computationally efficient and reasonably accurate BLPnet model for a new ALPR system.
 The model is designed to efficiently and correctly return the VLP number by considering a three\hyp phase top\hyp to\hyp bottom approach.
 This means, the model starts from detecting the vehicle first, and then any possible VLP followed by the recognition of characters in the VLP.
This consideration is observed to be more  effective with lower-computational time and accuracy than the prominent ALPR systems.\looseness -1

Moreover,  the model generates the actual license number of the vehicle from the recognized characters using our simple, yet effective mapping algorithm with a set of predefined cases of registration area\hyp codes.
The model also performed well without compromising accuracy to tackle challenging conditions that cause rotated,  blurry or noisy frames at the input.
The preliminary results with the reasonably faster and accurate performance of the model suggest that it would be promising for the real\hyp time ALPR application with a continuing development in future.\looseness -1

 \bibliographystyle{elsarticle-num}


\section{introduction}
Identifying infected people is a critical step in dealing with pandemics caused by viral diseases. However, testing a large number of people individually might be prohibitively expensive for practical reasons. For this reason, we need to deploy strategies that allow efficient testing.
 \emph{Group Testing} (GT) has been shown as an efficient strategy in reducing the number of tests required to test for pandemics.
An instance of the GT problem includes a  set $\mbf{S}$ of $N$ individuals which includes a small subset of infected individuals. The GT procedure consists of a sequence of tests, such that each test indicates whether there are one or more infected individuals in a given subset of $\mbf{S}$. The goal of the GT procedure is to identify the subset of infected individuals through the minimum number of tests.

The GT problem has been the subject of many studies. Most studies have focused on the following two models \cite{8989403,8635972,DBLP:journals/corr/abs-2110-10110,MCDERMOTT2021532}: (i) a  \emph{combinatorial} model which assumes that the number of infected individuals is fixed and known; (ii) a \emph{probabilistic} model which assumes that each individual is infected with a certain probability. There are also two types of GT algorithms:  \emph{non-adaptive}, and \emph{adaptive}.
In this paper, we are interested in {non-adaptive} GT strategies, where all tests are designed in advance. This is in contrast to {adaptive} strategies, in which the design of each test depends on the results of the previous tests~\cite{9646274,8437774,8919896,CIT-099}.

Motivated by practical scenarios where the outcome of the tests can be affected by noise, we focus on the \emph{noisy GT} setting, in which the outcome of a test can be flipped with some probability. In the noisy GT setting, the goal is to identify the set of infected individuals with high probability $(1-\varepsilon)$, for small values of $\varepsilon$. We also focus on a variation of a probabilistic GT model in which the prior infection probability is not uniform and in which there is a certain amount of side information on the distribution of the infected individuals available to the GT algorithm.

A GT algorithm consists of two parts: \emph{encoding} and \emph{decoding}. The encoding part is concerned with the test design, i.e., the decision on which individuals to include in each test. The decoding part is concerned with identifying the infected individuals given the test design and outcomes of the tests. Different decoding algorithms such as linear programming, combinatorial orthogonal matching pursuit, definite defectives, belief propagation (BP), and separate decoding of items have been proposed for noisy non-adaptive GT. A thorough review and comparison of these algorithms is provided in~\cite{CIT-099}.
In the context of GT, it is extremely difficult to analyze the performance of BP algorithms even for the asymptotic regime. To the best of our knowledge, no theoretical analysis has been provided for the BP-based GT algorithms so far.
However, empirical evidence suggests that the BP algorithm results in lower error probabilities compared to other algorithms for the probabilistic model.
BP is a message passing algorithm that passes messages over the edges in the underlying factor graph representation of the GT problem. For a cycle-free factor graph, BP decoding is equivalent to Maximum a Posteriori (MAP) decoding. However, in the presence of loops in the factor graph, BP becomes suboptimal.

This paper focuses on leveraging the side information for improving the performance of BP-based decoding algorithms for noisy GT. In the context of testing for viral infections, different forms of side information can be exploited including the prevalence rate, individuals' symptoms, family structure, community structure, and contact tracing information. It has been shown that side information can be used to reduce the required tests~\cite{zhu2020noisy,DBLP:journals/corr/abs-2007-08111,DBLP:journals/corr/abs-2012-02804,DBLP:journals/corr/abs-2106-02699,DBLP:journals/corr/abs-2101-02405}.
   For example, Zhu \emph{et al.}~\cite{zhu2020noisy} show that the number of tests can be reduced if the prior information about the prevalence rate is takend into account. Nikolopoulos \emph{et al.}~\cite{DBLP:journals/corr/abs-2007-08111,DBLP:journals/corr/abs-2012-02804} and Ahn \emph{et al.}~\cite{DBLP:journals/corr/abs-2101-02405} show that utilizing community structure also leads to a lower number of tests.
While the focus of these works is on the encoder design, in our work we focus on leveraging the side information for the efficient decoder design.

\subsection{Contribution.}

In this work, first, we propose a probabilistic model, referred to as an \emph{interaction model}, that captures the side information about the probability distribution of the infected individuals. Our model is motivated by the availability of contact tracing information which can be collected from surveys and mobile phone applications~\cite{KleinmanE653,doi:10.1080/19371918.2020.1806170,munzert2021tracking}.
 Next, we present a decoding scheme, based on belief propagation, that leverages the interaction model to improve the decoding accuracy. Our results indicate that the proposed algorithm achieves higher success probability and lower false-negative and false-positive rates when compared to the traditional belief propagation especially in the high noise regime.

\section{Probabilistic Model}\label{sec:SN}
Throughout the paper, we denote vectors and matrices by bold-face small and capital letters, respectively. For an integer $i\geq 1$, we denote $\{1,\dots,i\}$ by $[i]$. Let $\mbf{S}$ be a set. The set of all subsets of size $\ell$ for set $\mbf{S}$ is denoted by $\{\mbf{S}\}_\ell$.

Our model assumes that there are two points in time, namely time $0$ and time $1$ such that time $0$ occurs prior to time $1$. Let $N$ be the total number of individuals.
We define the vector $\mathbf{x}^{(0)}\in \{0,1\}^N$ to represent the status of $N$ individuals at time $0$, such that $x_i^{(0)}$ is $1$ if the $i$-th individual is infected at time $0$, and is $0$ otherwise.
Similarly, we define the vector $\mathbf{x}^{(1)}\in \{0,1\}^N$ to represent the status of $N$ individuals at time $1$. We assume that at time 0, the probability of an individual being infected is equal to the prevalence rate $p$, and that the infection of each individual at time $0$ occurs independently of other individuals. The probability of the individual to be infected at time 1 depends on their probability to be infected at time 0 as well as their interaction with other individuals. The interactions of individuals between time 0 and time 1 is captured by the \emph{interaction graph}. For each individual $i$, the graph includes nodes $x_i^{(0)}$ and  $x_i^{(1)}$ that represent that individual at times 0 and 1, respectively. For each individual $i$, the graph has an \emph{interaction node} $I_i$  that captures interactions between individual $i$ and other individuals from time $0$ to time $1$.  In particular, the graph contains an edge $(x_j^{(0)},I_i)$  for each individual $j$ who have been in contact with individual $i$ from time $0$ to time $1$. An example of an interaction graph is shown in Fig.~\ref{fig:Interaction}.

\begin{figure}
\centering
\begin{tikzpicture}
\def\horzgap{2in};
\def \gapVN{0.45in};
\def \gapCN{0.45in};

\def\nodewidth{0.25in};
\def\nodewidthA{0.3in};
\def \edgewidth{0.01in};
\def\ext{0.1in};
\def \dotwidth{0.3mm};

\tikzstyle{dots} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=black, minimum height=\dotwidth, minimum width=\dotwidth]
\tikzstyle{check} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=white, minimum height=0.75*\nodewidthA, minimum width=0.75*\nodewidthA]
\tikzstyle{check1} = [rectangle, draw, dotted,line width=0.4mm,  inner sep=0mm, fill=white, minimum height=0.6*\nodewidthA, minimum width=0.6*\nodewidthA]
\tikzstyle{bit0} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=white, minimum size=\nodewidth]
\tikzstyle{bit1} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=blue, minimum size=\nodewidth]
\tikzstyle{bituncover} = [circle, draw=none, line width=0.2mm, inner sep=0mm, fill=gray, minimum size=\nodewidthA]
\tikzstyle{edgesock} = [circle, inner sep=0mm, minimum size=\edgewidth,draw, fill=white]

\foreach \vn in {1,2,3,4,5}{
 \node[bit0] (vn0\vn) at (-0.7*\horzgap,-\vn*\gapVN) {};
}
\path (vn01) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{1}^{(0)}$}};
\path (vn02) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{2}^{(0)}$}};
\path (vn03) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{3}^{(0)}$}};
\path (vn04) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{4}^{(0)}$}};
\path (vn05) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{5}^{(0)}$}};

\foreach \vn in {1,2,3,4,5}{
 \node[bit0] (vn1\vn) at (0.5*\horzgap,-\vn*\gapVN) {};
}
\path (vn11) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{1}^{(1)}$}};
\path (vn12) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{2}^{(1)}$}};
\path (vn13) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{3}^{(1)}$}};
\path (vn14) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{4}^{(1)}$}};
\path (vn15) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{5}^{(1)}$}};

\foreach \cn in {1,...,5}{
\node[check] (cn\cn) at (0,-\cn*\gapVN) {};
}
\draw[line width=0.2mm] (vn01.east)--(cn1.west);
\draw[line width=0.2mm] (vn03.east)--(cn1.west);
\draw[line width=0.2mm] (vn05.east)--(cn1.west);

\draw[line width=0.2mm] (vn02.east)--(cn2.west);

\draw[line width=0.2mm] (vn01.east)--(cn3.west);
\draw[line width=0.2mm] (vn03.east)--(cn3.west);

\draw[line width=0.2mm] (vn04.east)--(cn4.west);
\draw[line width=0.2mm] (vn05.east)--(cn4.west);

\draw[line width=0.2mm] (vn01.east)--(cn5.west);
\draw[line width=0.2mm] (vn04.east)--(cn5.west);
\draw[line width=0.2mm] (vn05.east)--(cn5.west);

\draw[line width=0.2mm] (cn1.east)--(vn11.west);
\draw[line width=0.2mm] (cn2.east)--(vn12.west);
\draw[line width=0.2mm] (cn3.east)--(vn13.west);
\draw[line width=0.2mm] (cn4.east)--(vn14.west);
\draw[line width=0.2mm] (cn5.east)--(vn15.west);

\def\moveX {3.2*\nodewidth};
\def\moveXA {2*\nodewidth};

\path (cn1)++(0,0) node ()[scale=0.25] {\Huge{$I_1$}};
    \path (cn2)++(0,0) node ()[scale=0.25] {\Huge{$I_2$}};
    \path (cn3)++(0,0) node ()[scale=0.25] {\Huge{$I_3$}};
    \path (cn4)++(0,0) node ()[scale=0.25] {\Huge{$I_4$}};
     \path (cn5)++(0,0) node ()[scale=0.25] {\Huge{$I_5$}};

\end{tikzpicture}
\caption{An example of an interaction graph. Nodes $x_i^{(0)}$ and  $x_i^{(1)}$ represent individual $i$ at times 0 and 1, respectively. An interaction node $I_i$ captures interactions between individual $i$ and other individuals from time $0$ to time $1$.} \label{fig:Interaction
\end{figure}

We assume that an infected individual infects a healthy individual with probability $q$, referred to as \emph{contagion probability.} It is also assumed that if an individual is infected at time $0$, they remain infected by time $1$. The interaction model can be used to find the probability that an individual at time $1$ is infected. The $i$th individual is not infected at time $1$ if the following holds: $1)$ the $i$th individual is not infected at time $0$, and $2)$ other individuals in contact with the $i$th individual either are not infected at time $0$ or, if infected, they do not infect the $i$th individual. Thus, the probability of individual $i$ to be infected at time 1 can be calculated as follows:
\[
P\Big(x_i^{(1)}=0\Big) =(1-p)\Big(1-p+p(1-q)\Big)^{d_i}
=(1-p)(1-pq)^{d_i},
 \]
where $d_i$ is the number of individuals which interact with individual $i$ from time $0$ to time $1$. The probability that individual $i$ is infected at time $1$ is given by
\begin{equation}\label{eq:prior}
\pi_i=P\Big(x_i^{(1)}=1\Big)=1-P\Big(x_i^{(1)}=0\Big) =1-(1-p)(1-pq)^{d_i}.
\end{equation}
Our ultimate goal is to test and identify infected individuals at time $1$, assuming the knowledge of the probabilistic model described above. This model can be extended to capture interactions between individuals in more than one round. For ease of exposition, we consider only one round of interactions in this study.

In non-adaptive group testing, designing a testing scheme consisting of $M$ tests is equivalent to the construction of a binary matrix with $M$ rows which is referred to as a testing matrix.
We let matrix ${\textbf{A}\in \{0,1\}^{M\times N}}$ denote the testing matrix. The entry $(t,i)$ of matrix ${\textbf{A}}$ is denoted by $a_{t,i}$.
If $a_{t,i}=1$, it means that the $i$-th item is present in the $t$-th test.
The design of testing matrices for group testing has been studied extensively (see e.g., \cite{CIT-099}).
Our proposed algorithms are applicable for any testing matrix; however, simulation results are presented for Bernoulli designs. In a Bernoulli design, each individual is included in each test independently at random with some fixed probability $\nu/K$ where $K$ is the average number of defective items and $\nu$ is a constant.

 The standard noiseless group testing is formulated component-wise using the Boolean OR operation as ${y_t=\bigvee_{i=1}^{N} a_{t,i}x_i^{(1)}}$ where $y_t$ and $\bigvee$ are the $t$th test result and a Boolean OR operation, respectively. In this paper, we consider the widely-adopted binary symmetric noise model where the values $\bigvee_{i=1}^{N} a_{t,i}x_i^{(1)}$ are flipped independently at random with a given probability. The $t$th test result in a binary symmetric noise model is given by
\[
     y_t=\begin{cases}
       \bigvee_{i=1}^{N} a_{t,i}x_i^{(1)} &\quad\text{with probability}~ 1-\rho,\\
        1\oplus \bigvee_{i=1}^{N} a_{t,i}x_i^{(1)} &\quad \text{with probability} ~\rho,\\

     \end{cases}
\]
where $\oplus$ is the XOR operation and $\rho$ is the probability that the values $\bigvee_{i=1}^{N} a_{t,i}x_i^{(1)}$ are flipped.
Note that this model and the proposed algorithms can be easily extended to include the general binary noise model where the values $\bigvee_{i=1}^{N} a_{t,i}x_i^{(1)}$ are flipped from $0$ to $1$ and from $1$ to $0$ with different probabilities independently.
However, for ease of exposition, we focus only on the binary symmetric noise model.
 We let vector $\mathbf{y}\in\{0,1\}^M$ denote the outcomes of the $M$ tests.

 Our objective is to design a decoding algorithm that performs well under the following three metrics: (i) success probability which captures the probability that all infected individuals are identified correctly; (ii) False-Negative Rate (FNR), which is the number of infected individuals falsely classified as healthy over the total number of infected individuals, (iii) False-Positive Rate (FPR), defined as the ratio of the number of healthy individuals falsely classified as infected and the total number of healthy individuals.

\section{Proposed Decoding Algorithms \label{sec:main results}}
In this section, we describe the proposed decoding algorithms for retrieving the status vector $\mathbf{x}^{(1)}$ from the test results vector $\mathbf{y}$ and the testing matrix ${\textbf{A}}$.
\subsection{Belief Propagation Using Initial Prior Probabilities}\label{subsec:bp}
Message passing algorithms are utilized to solve inference problems, optimization problems, and constraint satisfaction problems. In an inference problem, there are some noisy measurements as input, and the goal is to infer the value of some unobserved variables from those measurements. It is impossible, in general, to make those inferences with complete certainty, but one can try to obtain the most probable value of the unobserved variables\cite{mackay2003information,russell2003artificial}. In a probabilistic noisy group testing, we intend to perform a Maximum a Posteriori (MAP) estimation to find the status vector $\hat{\mbf{x}}^{(1)}$ given the test results vector $\mbf{u}$.
\begin{equation}\label{eq:cost}
\mathop{\mathrm{arg\,max}_{\hat{\mbf{x}}^{(1)}} P\Big(\mbf{x}^{(1)}=\hat{\mbf{x}}^{(1)}\Big)P\Big(\mbf{y}=\mbf{u}\vert \mbf{x}^{(1)}=\hat{\mbf{x}}^{(1)}\Big)
\end{equation}
This problem can be solved using exhaustive search when the vector $\mbf{x}^{(1)}$ is small, e.g., $N\approx 10-20$.
However, the exhaustive search approach rapidly becomes intractable when $N$ increases. For instance, when $N=100$ (a relatively small problem), the number of different configurations for the vector $\mbf{x}^{(1)}$ is $2^{100}$.

An alternative solution is to find the marginals of the posterior distribution for each item using Belief Propagation (BP). BP is a message passing algorithm for performing inference on factor graphs with the purpose of calculating the marginal distribution for each unobserved variable, conditional on observed variables. A factor graph is a type of probabilistic graphical model which is used to visualize and precisely define the underlying optimization problem. Factor graphs are bipartite graphs with two types of nodes referred to as variable nodes and factor nodes. See Fig.~\ref{fig:Bipartite} for an example of a factor graph.  The variable nodes which represent the variables in the optimization problem are represented by circles. The factor nodes show how the overall cost function can be factorized into local cost functions and are represented by squares. There is an edge between the variables that are involved in a local cost function and the factor node representing that local cost function.

We assume that the only side information we have is the prevalence rate at time $0$. Since we have no information about the status of individuals at time $1$, we consider the prior probability of each individual being infected at time $1$ to be equal to the prevalence rate $p$, independent of other individuals. The overall cost function in~\eqref{eq:cost} can be factorized as follows.
\begin{equation}
 P\Big(\mbf{x}^{(1)}=\hat{\mbf{x}}^{(1)}\Big)P\Big(\mbf{y}=\mbf{u}\vert \mbf{x}^{(1)}=\hat{\mbf{x}}^{(1)}\Big) =\Bigg[\prod_{i=1}^N P\Big({x}_i^{(1)}=\hat{{x}}_i^{(1)}\Big)\Bigg] \times \Bigg[\prod_{t=1}^M P\Big(y_t=u_t \Big \vert \Big\{{{x}}_{i}^{(1)}=\hat{{x}}_{i}^{(1)}\Big\}_{i\in \mathcal{N}(t)}\Big)\Bigg],
\end{equation}
where $\mathcal{N}(t)$ denotes the indices of individuals involved in the $t$th test.

In order to apply BP, we consider the factor graph (Tanner graph) representation of the group testing scheme. In the Tanner graph, there are $N$ variable nodes that represent individuals at time $1$. There are also $M$ factor nodes that represent the tests. Each test factor node corresponds to the conditional probability distribution of a test result, given the observed variable nodes.
Each individual in the Tanner graph is connected to the test in which the individual participates, according to the testing matrix. There are also $N$ factor nodes that correspond to the \textit{a priori} probability distribution of the variable nodes. Since these factor nodes are usually not exhibited in a Tanner graph, we show them using dotted squares in Fig.~\ref{fig:Bipartite}.

\begin{figure}
\centering
\begin{tikzpicture}
\def\horzgap{2in};
\def \gapVN{0.45in};
\def \gapCN{0.45in};

\def\nodewidth{0.25in};
\def\nodewidthA{0.3in};
\def \edgewidth{0.01in};
\def\ext{0.1in};
\def \dotwidth{0.3mm};

\tikzstyle{dots} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=black, minimum height=\dotwidth, minimum width=\dotwidth]
\tikzstyle{check} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=white, minimum height=0.75*\nodewidthA, minimum width=0.75*\nodewidthA]
\tikzstyle{check1} = [rectangle, draw, dotted,line width=0.4mm,  inner sep=0mm, fill=white, minimum height=0.6*\nodewidthA, minimum width=0.6*\nodewidthA]
\tikzstyle{bit0} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=white, minimum size=\nodewidth]
\tikzstyle{bit1} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=blue, minimum size=\nodewidth]
\tikzstyle{bituncover} = [circle, draw=none, line width=0.2mm, inner sep=0mm, fill=gray, minimum size=\nodewidthA]
\tikzstyle{edgesock} = [circle, inner sep=0mm, minimum size=\edgewidth,draw, fill=white]

\foreach \vn in {1,2,3,4,5}{
 \node[bit0] (vn\vn) at (0,-\vn*\gapVN) {};
}
\path (vn1) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{1}^{(1)}$}};
\path (vn2) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{2}^{(1)}$}};
\path (vn3) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{3}^{(1)}$}};
\path (vn4) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{4}^{(1)}$}};
\path (vn5) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{5}^{(1)}$}};
\foreach \cn in {1,...,3}{
\node[check] (cn\cn) at (0.7*\horzgap,-\cn*\gapCN-0.4in) {};
}
\draw[line width=0.2mm] (vn1.east)--(cn1.west);
\draw[line width=0.2mm] (vn3.east)--(cn1.west);
\draw[line width=0.2mm] (vn5.east)--(cn1.west);
\draw[line width=0.2mm] (vn2.east)--(cn2.west);
\draw[line width=0.2mm] (vn3.east)--(cn2.west);
\draw[line width=0.2mm] (vn4.east)--(cn2.west);
\draw[line width=0.2mm] (vn2.east)--(cn3.west);
\draw[line width=0.2mm] (vn4.east)--(cn3.west);
\draw[line width=0.2mm] (vn5.east)--(cn3.west);
\foreach \pn in {1,...,5}{\node[check1] (pn\pn) at (-0.7in,-\pn*\gapVN) {};}
\draw[dotted, line width=0.4mm] (pn1.east)--(vn1.west);
\draw[dotted, line width=0.4mm] (pn2.east)--(vn2.west);
\draw[dotted, line width=0.4mm] (pn3.east)--(vn3.west);
\draw[dotted, line width=0.4mm] (pn4.east)--(vn4.west);
\draw[dotted, line width=0.4mm] (pn5.east)--(vn5.west);
\def\moveX {3.2*\nodewidth};
\def\moveXA {2*\nodewidth};
\path (cn1)++(0,0) node ()[scale=0.25] {\Huge{$y_1$}};
    \path (cn2)++(0,0) node ()[scale=0.25] {\Huge{$y_2$}};
    \path (cn3)++(0,0) node ()[scale=0.25] {\Huge{$y_3$}};
\path (pn1.north)++(0,0.5) node ()[scale=0.3] {\Huge{Priors}};
\path (vn1.north)++(0,0.5) node ()[scale=0.3] {\Huge{Individuals}};
\path (cn1.north)++(0,0.5) node ()[scale=0.3] {\Huge{Tests}};
\end{tikzpicture}
\caption{An example of a factor graph representing a group testing scheme. A variable node $x_{i}^{(1)}$ represents individual $i$ at time $1$. A factor node $y_t$ represents test $t$. Factor nodes represented by dotted squares correspond to the \textit{a priori} probability distribution of the variable nodes. } \label{fig:Bipartite
\end{figure}

For a cycle-free factor graph, BP decoding is equivalent to Maximum a Posteriori (MAP) decoding. However, in the presence of loops in the factor graph, BP becomes suboptimal. In other words, loopy BP computes an approximation of the marginals of the posterior distribution for each variable node. For a loopy BP algorithm, the messages are passed iteratively from variable nodes to factor nodes and vice versa. We let ${\mu}_{i \rightarrow t} = [\mu_{i\rightarrow t}(0) \ \mu_{i\rightarrow t}(1)]$ and ${\mu}_{t\rightarrow i} = [\mu_{t \rightarrow i}(0) \ \mu_{t\rightarrow i}(1)]$ denote the message from individual $i$ to test $t$ and the message from test $t$ to individual $i$, respectively. In general, the message from variable node $i$ to factor node $t$ is given by computing the product of all incoming messages from the neighboring factor nodes of variable node $i$ excluding the message from factor node $t$.
\begin{equation}\label{eq:itm-to-tst}
\begin{cases}
 \mu_{i\rightarrow t}(0) \propto (1-p) \displaystyle \prod_{t^{\prime}\in \mathcal{N}(i)\setminus \{t\}}\mu_{t^{\prime}\rightarrow i}(0), \\
 \mu_{i\rightarrow t}(1) \propto p \displaystyle \prod_{t^{\prime}\in \mathcal{N}(i)\setminus \{t\}}\mu_{t^{\prime}\rightarrow i}(1),\\
  \end{cases}
\end{equation}
where $\propto$ indicates equality up to a
normalizing constant, and $\mathcal{N}(i)$ denotes the indices of tests in which item $i$ participates. Note that these messages are probability distributions on \{0,1\}, i.e., $\mu_{i\rightarrow t}(0) + \mu_{i\rightarrow t}(1)=1 $.
Since we assume that the prior probability of each individual being infected at time $1$ is equal to $p$, the messages are initialized by
\begin{equation}\label{eq:intl}
 \mu_{i\rightarrow t}(1)=1-\mu_{i\rightarrow t}(0)=p.
\end{equation}
The messages from factor nodes to variable nodes are computed as follows. The message from test $t$ to individual $i$ is given by
\begin{equation}\label{eq:ctov}
 \mu_{t\rightarrow i}(\hat{{x}}_{i}^{(1)})=  \sum_{\substack{\hat{{x}}_{i^\prime}^{(1)}\in\{0,1\},~{i^{\prime}\in \mathcal{N}(t)\setminus\{i\}}\\u_t\in \{0,1\}}} \Bigg[P\Big(y_t=u_t\Big\vert {x}_{i}^{(1)}=\hat{{x}}_{i}^{(1)},\Big\{{x}_{i^\prime}^{(1)}=\hat{{x}}_{i^\prime}^{(1)}\Big\}_{i^{\prime}\in \mathcal{N}(t)\setminus\{i\}}\Big) \prod_{i^{\prime}\in \mathcal{N}(t)\setminus\{i\}}\mu_{i^{\prime}\rightarrow t}\Big(\hat{{x}}_{i^\prime}^{(1)}\Big)\Bigg].
\end{equation}
As was shown in~\cite{barbier2020strong,2018,coja2018belief,doi:10.1080/00018732.2016.1211393}, the equation~\eqref{eq:ctov} can be simplified as follows. If $y_t=0$, we have
\begin{equation}\label{eq:tst-to-itm0}
\begin{cases}
 \mu_{t\rightarrow i}(0) \propto  \rho + (1-2\rho) \displaystyle \prod_{i^{\prime}\in \mathcal{N}(t)\setminus \{i\}}\mu_{i^{\prime}\rightarrow t}(0), \\
 \mu_{t\rightarrow i}(1) \propto  \rho, \\
  \end{cases}
\end{equation} and if $y_t=1$, we have
\begin{equation}\label{eq:tst-to-itm1}
\begin{cases}
 \mu_{t\rightarrow i}(0) \propto  1-\rho - (1-2\rho) \displaystyle \prod_{i^{\prime}\in \mathcal{N}(t)\setminus \{i\}}\mu_{i^{\prime}\rightarrow t}(0), \\
 \mu_{t\rightarrow i}(1) \propto  1-\rho. \\
  \end{cases}
\end{equation}

We perform a fixed point iteration using the BP equations~\eqref{eq:itm-to-tst},~\eqref{eq:tst-to-itm0}, and~\eqref{eq:tst-to-itm1}. We stop the algorithm after a fixed number $T$ of iterations. The parameter $T$ is chosen experimentally. In the end, we compute the marginals of the posterior distribution for each variable node by computing the product of all incoming messages from the neighboring factor nodes of that variable node.

\begin{equation}
\begin{cases}
{{\phi}}\Big(x_{i}^{(1)}=0\Big) \propto  (1-p) \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(0), \\
  {{\phi}}\Big(x_{i}^{(1)}=1\Big) \propto  p \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(1). \\
  \end{cases} \vspace{0.1cm}
\end{equation}
For convenience, we work with the Log-Likelihood Ratio (LLR) of a marginal defined as
\begin{equation}\label{eq:llr}
\lambda_i=\ln \frac{ {{\phi}}\Big(x_{i}^{(1)}=1\Big)}{{{\phi}}\Big(x_{i}^{(1)}=0\Big)}=\ln \frac{p}{1-p} + \sum_{t\in \mathcal{N}(i)}\ln \frac{\mu_{t\rightarrow i}(1)}{\mu_{t\rightarrow i}(0)}.
\end{equation}

We consider a threshold, $\tau$, and announce the $i$th individual infected if $ \lambda_i \geq \tau$. A natural threshold one can choose is $\tau=0$. Note that values other than $0$ are also permissible.
Algorithm~\ref{alg:BP} defines the belief propagation using initial prior probabilities algorithm.

\begin{algorithm}
\caption{Belief Propagation Using Initial Prior Probabilities}
\label{alg:BP}
\label{CHalgorithm}
\begin{algorithmic}[1]
\State Initialize $\mu_{i\rightarrow t}(1)=1-\mu_{i\rightarrow t}(0)=p$ $\forall i\in [N], \forall t\in \mathcal{N}(i)$
\For{$\ell=1,2,\cdots,T$ }
\State Compute $\mu_{t\rightarrow i}(0)$ and $\mu_{t\rightarrow i}(1)$ $\forall t\in [M], \forall i\in \mathcal{N}(t)$ using~\eqref{eq:tst-to-itm0} and~\eqref{eq:tst-to-itm1}
\State Compute $\mu_{i\rightarrow t}(0)$ and $\mu_{i\rightarrow t}(1)$ $\forall i\in [N], \forall t\in \mathcal{N}(i)$ using~\eqref{eq:itm-to-tst}
\EndFor
\State Compute $\lambda_i$ $\forall  i\in [N]$ using~\eqref{eq:llr}
\end{algorithmic}
\end{algorithm}

\subsection{Belief Propagation Using Updated Prior Probabilities.}

In this scheme, instead of using the prevalence rate at time $0$ for the probability that an individual is infected at time $1$, we use the updated prior probability $\pi_i$, $i\in [N]$, given by~\eqref{eq:prior}. We perform the BP algorithm in Section~\ref{subsec:bp} where in the equations~\eqref{eq:itm-to-tst}-\eqref{eq:llr}, the initial prior probability $p$ is replaced by the updated prior probability $\pi_i$, for each $i\in [N]$.

\begin{figure}
\centering
\begin{tikzpicture}
\def\horzgap{2in};
\def \gapVN{0.45in};
\def \gapCN{0.45in};

\def\nodewidth{0.25in};
\def\nodewidthA{0.3in};
\def \edgewidth{0.01in};
\def\ext{0.1in};
\def \dotwidth{0.3mm};

\tikzstyle{dots} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=black, minimum height=\dotwidth, minimum width=\dotwidth]
\tikzstyle{check} = [rectangle, draw,line width=0.2mm,  inner sep=0mm, fill=white, minimum height=0.75*\nodewidthA, minimum width=0.75*\nodewidthA]
\tikzstyle{check1} = [rectangle, draw, dotted,line width=0.4mm,  inner sep=0mm, fill=white, minimum height=0.6*\nodewidthA, minimum width=0.6*\nodewidthA]
\tikzstyle{bit0} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=white, minimum size=\nodewidth]
\tikzstyle{bit1} = [circle, draw, line width=0.2mm, inner sep=0mm, fill=blue, minimum size=\nodewidth]
\tikzstyle{bituncover} = [circle, draw=none, line width=0.2mm, inner sep=0mm, fill=gray, minimum size=\nodewidthA]

\tikzstyle{edgesock} = [circle, inner sep=0mm, minimum size=\edgewidth,draw, fill=white]

\foreach \vn in {1,2,3,4,5}{
 \node[bit0] (vn0\vn) at (-0.65*\horzgap,-\vn*\gapVN) {};
}
\path (vn01) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{1}^{(0)}$}};
\path (vn02) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{2}^{(0)}$}};
\path (vn03) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{3}^{(0)}$}};
\path (vn04) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{4}^{(0)}$}};
\path (vn05) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{5}^{(0)}$}};

\foreach \vn in {1,2,3,4,5}{
 \node[bit0] (vn1\vn) at (0.1*\horzgap,-\vn*\gapVN) {};
}
\path (vn11) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{1}^{(1)}$}};
\path (vn12) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{2}^{(1)}$}};
\path (vn13) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{3}^{(1)}$}};
\path (vn14) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{4}^{(1)}$}};
\path (vn15) ++(0,0) node()[scale=0.25, inner sep=0mm] {\Huge{$x_{5}^{(1)}$}};

\foreach \cn in {1,...,5}{
\node[check] (cn\cn) at (-0.15*\horzgap,-\cn*\gapVN) {};
}
\draw[line width=0.2mm] (vn01.east)--(cn1.west);
\draw[line width=0.2mm] (vn03.east)--(cn1.west);
\draw[line width=0.2mm] (vn05.east)--(cn1.west);

\draw[line width=0.2mm] (vn02.east)--(cn2.west);

\draw[line width=0.2mm] (vn01.east)--(cn3.west);
\draw[line width=0.2mm] (vn03.east)--(cn3.west);

\draw[line width=0.2mm] (vn04.east)--(cn4.west);
\draw[line width=0.2mm] (vn05.east)--(cn4.west);

\draw[line width=0.2mm] (vn01.east)--(cn5.west);
\draw[line width=0.2mm] (vn04.east)--(cn5.west);
\draw[line width=0.2mm] (vn05.east)--(cn5.west);

\draw[line width=0.2mm] (cn1.east)--(vn11.west);
\draw[line width=0.2mm] (cn2.east)--(vn12.west);
\draw[line width=0.2mm] (cn3.east)--(vn13.west);
\draw[line width=0.2mm] (cn4.east)--(vn14.west);
\draw[line width=0.2mm] (cn5.east)--(vn15.west);

\foreach \cn in {1,...,3}{
\node[check] (cn1\cn) at (0.6*\horzgap,-\cn*\gapCN-0.4in) {};
}

\draw[line width=0.2mm] (vn11.east)--(cn11.west);
\draw[line width=0.2mm] (vn13.east)--(cn11.west);
\draw[line width=0.2mm] (vn15.east)--(cn11.west);

\draw[line width=0.2mm] (vn13.east)--(cn12.west);
\draw[line width=0.2mm] (vn14.east)--(cn12.west);
\draw[line width=0.2mm] (vn11.east)--(cn12.west);

\draw[line width=0.2mm] (vn12.east)--(cn13.west);
\draw[line width=0.2mm] (vn14.east)--(cn13.west);
\draw[line width=0.2mm] (vn15.east)--(cn13.west);
\draw[line width=0.2mm] (vn11.east)--(cn13.west);

\foreach \pn in {1,...,5}{\node[check1] (pn\pn) at (-0.9*\horzgap,-\pn*\gapVN) {};}
\draw[dotted, line width=0.4mm] (pn1.east)--(vn01.west);
\draw[dotted, line width=0.4mm] (pn2.east)--(vn02.west);
\draw[dotted, line width=0.4mm] (pn3.east)--(vn03.west);
\draw[dotted, line width=0.4mm] (pn4.east)--(vn04.west);
\draw[dotted, line width=0.4mm] (pn5.east)--(vn05.west);

\def\moveX {3.2*\nodewidth};
\def\moveXA {2*\nodewidth};

\path (cn1)++(0,0) node ()[scale=0.25] {\Huge{$I_1$}};
    \path (cn2)++(0,0) node ()[scale=0.25] {\Huge{$I_2$}};
    \path (cn3)++(0,0) node ()[scale=0.25] {\Huge{$I_3$}};
    \path (cn4)++(0,0) node ()[scale=0.25] {\Huge{$I_4$}};
     \path (cn5)++(0,0) node ()[scale=0.25] {\Huge{$I_5$}};

   \path (cn11)++(0,0) node ()[scale=0.25] {\Huge{$y_1$}};
    \path (cn12)++(0,0) node ()[scale=0.25] {\Huge{$y_2$}};
    \path (cn13)++(0,0) node ()[scale=0.25] {\Huge{$y_3$}};

\end{tikzpicture}

\caption{An example of a combined graph. Nodes $x_i^{(0)}$ and  $x_i^{(1)}$ represent individual $i$ at times 0 and 1, respectively. An interaction node $I_i$ captures interactions between individual $i$ and other individuals from time $0$ to time $1$. A factor node $y_t$ represents test $t$. Factor nodes represented by dotted squares correspond to the \textit{a priori} probability distribution of the status of individuals at time $0$.} \label{fig:Comb
\end{figure}

 \subsection{Belief Propagation on Combined Graphs}
In this scheme, assuming that the contact tracing information is available, we form the interaction graph and combine it with the Tanner graph corresponding to the testing matrix. An example of a combined graph is presented in Fig.~\ref{fig:Comb}.
We then perform a BP algorithm over the combined graph. Note that there are two sets of variable nodes in the combined graph, $\Big\{x_i^{(0)}\Big\}_{i\in[N]}$ and $\Big\{x_i^{(1)}\Big\}_{i\in[N]}$. We are interested in computing the marginals of the posterior distribution for $\Big\{x_i^{(1)}\Big\}_{i\in[N]}$. There are also three different types of factor nodes. The interaction node $I_i$ corresponds to the conditional probability that individual $i$ at time $1$ is infected or not, given the status of individuals at time $0$ who have been in contact with individual $i$. The test factor node $y_t$ corresponds to the conditional probability that the result of test $t$ is equal to a one or zero, given the status
of neighboring individuals at time $1$. Furthermore, there are $N$ factor nodes, represented by dashed squares, that correspond to the \textit{a priori} probability that each of the individuals at time $0$ is infected or not. The combined graph in fact represents the factorization in the following joint probability mass function.

\begin{multline*}
P\Big(\mbf{x}^{(0)}=\hat{\mbf{x}}^{(0)},\mbf{x}^{(1)}=\hat{\mbf{x}}^{(1)},\mbf{y}=\mbf{u}\Big)\\
=\Bigg[\prod_{i=1}^N P\Big({x}_i^{(0)}=\hat{{x}}_i^{(0)}\Big)\Bigg]\times \Bigg[\prod_{i=1}^N P\Big({x}_i^{(1)}=\hat{{x}}_i^{(1)}\Big \vert \Big\{{x}_{i^\prime}^{(0)}=\hat{{x}}_{i^\prime}^{(0)}\Big\}_{i^{\prime}\in \mathcal{N}(I_i)}\Big)\Bigg]\times \Bigg[\prod_{t=1}^M P\Big(y_t=u_t\Big \vert \Big \{{x}_{i}^{(1)}=\hat{{x}}_{i}^{(1)}\Big\}_{i\in \mathcal{N}(t)}\Big)\Bigg],
\end{multline*}
where $\mathcal{N}(I_i)$ denotes the indices of individuals at time $0$ who are connected to interaction node $I_i$.

As it has been mentioned before, in a loopy BP algorithm, the messages are passed iteratively from variable nodes to factor nodes and vice versa. In what follows, we show the flow of messages in one iteration. First, individuals at time $0$ send their messages to interaction nodes. We let ${\gamma}_{i \rightarrow I_j} = [\gamma_{i\rightarrow I_j}(0) \ \gamma_{i\rightarrow I_j}(1)]$ and ${\gamma}_{I_j \rightarrow i} = [\gamma_{I_j\rightarrow i}(0) \ \gamma_{I_j\rightarrow i}(1)]$ denote the message from individual $i$ at time $0$ to interaction node $I_j$ and the message from interaction node $I_j$ to individual $i$ at time $0$, respectively. It is easy to show that the messages ${\gamma}_{i \rightarrow I_j}$ can be computed as follows:
\begin{equation}\label{eq:itm0-to-fctr}
\begin{cases}
 \gamma_{i\rightarrow I_j}(0) \propto (1-p) \displaystyle \prod_{j^{\prime}\in \mathcal{N}(i)\setminus \{j\}}\gamma_{I_{j^{\prime}}\rightarrow i}(0), \\
 \gamma_{i\rightarrow I_j}(1) \propto p \displaystyle \prod_{j^{\prime}\in \mathcal{N}(i)\setminus \{j\}}\gamma_{I_{j^\prime}\rightarrow i}(1),\\
  \end{cases}
\end{equation}
where these messages are initialized by $\gamma_{i\rightarrow I_j}(1)=1-\gamma_{i\rightarrow I_j}(0)=p$. Then, the interaction nodes send their messages to individuals at time $1$.
We denote the message from interaction node $I_j$ to individual $j$ at time $1$ and the message from individual $j$ at time $1$ to interaction node $I_j$ by ${\delta}_{I_j \rightarrow j} = [\delta_{I_j\rightarrow j}(0) \ \delta_{I_j\rightarrow j}(1)]$ and ${\delta}_{j \rightarrow I_j} = [\delta_{j\rightarrow I_j}(0) \ \delta_{j\rightarrow I_j}(1)]$, respectively. It can be shown that the message ${\delta}_{I_j \rightarrow j}$ is given by

\begin{equation}\label{eq:Intr-to-itm1}
\begin{cases}

 \delta_{I_j\rightarrow j}(0)  \propto  \gamma_{j\rightarrow I_j}(0)  \displaystyle \prod_{i\in \mathcal{N}(I_j)\setminus \{j\}}\Big(1-q\gamma_{i\rightarrow I_j}(1)\Big ),\\

 \delta_{I_j\rightarrow j}(1)  \propto \gamma_{j\rightarrow I_j}(1)
  - \gamma_{j\rightarrow I_j}(0) \displaystyle \sum_{\ell=1}^{\vert \mathcal{N}(I_j) \vert -1}\sum_{S\in \Big\{\mathcal{N}(I_j)\setminus \{j\}\Big\}_\ell}\Big(-q\Big)^\ell\prod_{i\in S}\gamma_{i\rightarrow I_j}(1).

  \end{cases}
\end{equation}

In the next step, individuals at time $1$ send their messages to test factor nodes. We let ${\mu}_{i \rightarrow t} = [\mu_{i\rightarrow t}(0) \ \mu_{i\rightarrow t}(1)]$ and ${\mu}_{t\rightarrow i} = [\mu_{t \rightarrow i}(0) \ \mu_{t\rightarrow i}(1)]$ denote the message from individual $i$ at time $1$ to test $t$ and the message from test $t$ to individual $i$ at time $1$, respectively.
\begin{equation}\label{eq:itm1-to-tst}
\begin{cases}
 \mu_{i\rightarrow t}(0) \propto \delta_{I_i\rightarrow i}(0) \displaystyle \prod_{t^{\prime}\in \mathcal{N}(i)\setminus \{t\}}\mu_{t^{\prime}\rightarrow i}(0), \\
 \mu_{i\rightarrow t}(1) \propto \delta_{I_i\rightarrow i}(1) \displaystyle \prod_{t^{\prime}\in \mathcal{N}(i)\setminus \{t\}}\mu_{t^{\prime}\rightarrow i}(1).\\
  \end{cases}
\end{equation}
Now, test nodes send their messages to individuals at time $1$. The message ${\mu}_{t \rightarrow i}$ is calculated in a similar way as in~\eqref{eq:tst-to-itm0} and~\eqref{eq:tst-to-itm1}. Next, individuals at time $1$ send their messages to interaction nodes. It can be shown that the message from individual $i$ at time $1$ to interaction node $I_i$ is given by
\begin{equation}\label{eq:itm1-to-fctr}
\begin{cases}
 \delta_{i\rightarrow I_i}(0) \propto \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(0), \\
 \delta_{i\rightarrow I_i}(1) \propto \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(1).\\
  \end{cases}
\end{equation}
 Finally, interaction nodes send their messages to individuals at time $0$. The message from interaction node $I_j$ to individual $j$ at time $0$ is given as follows.
\begin{equation}\label{eq:fctr-to-itm0j}
\begin{cases}
 \gamma_{I_j\rightarrow j}(0) \propto \delta_{j\rightarrow I_j}(0) \displaystyle \prod_{i\in \mathcal{N}(I_j)\setminus \{j\}}\Big(1-q\gamma_{i\rightarrow I_j}(1)\Big )
 - \delta_{j\rightarrow I_j}(1) \displaystyle \sum_{\ell=1}^{\vert \mathcal{N}(I_j) \vert -1}\sum_{S\in \Big\{\mathcal{N}(I_j)\setminus \{j\}\Big\}_\ell}\Big(-q\Big)^\ell\prod_{i\in S}\gamma_{i\rightarrow I_j}(1), \\
 \gamma_{I_j\rightarrow j}(1) \propto \delta_{j\rightarrow I_j}(1).\\
  \end{cases}
\end{equation}
The message from interaction node $I_j$ to individual $i$ at time $0$, where $i\in \mathcal{N}(I_j)\setminus \{j\}$, is given by
\begin{equation}\label{eq:fctr-to-itm0i}
\begin{cases}
\begin{split}
 \gamma_{I_j\rightarrow i}(0) & \propto \gamma_{j\rightarrow I_j}(1) \delta_{j\rightarrow I_j}(1) + \gamma_{j\rightarrow I_j}(0) \delta_{j\rightarrow I_j}(0) \displaystyle \prod_{i^{\prime}\in \mathcal{N}(I_j)\setminus \{i,j\}}\Big(1-q\gamma_{i^{\prime}\rightarrow I_j}(1)\Big )\\
 & - \gamma_{j\rightarrow I_j}(0) \delta_{j\rightarrow I_j}(1)\sum_{\ell=1}^{\vert \mathcal{N}(I_j) \vert -2}\sum_{S\in \Big\{\mathcal{N}(I_j)\setminus \{i,j\}\Big\}_\ell}\Big(-q\Big)^\ell\prod_{i^{\prime}\in S}\gamma_{i^{\prime}\rightarrow I_j}(1),
\end{split}
 \\
\begin{split}
 \gamma_{I_j\rightarrow i}(1) & \propto \gamma_{j\rightarrow I_j}(1) \delta_{j\rightarrow I_j}(1) + (1-q)\gamma_{j\rightarrow I_j}(0) \delta_{j\rightarrow I_j}(0) \displaystyle \prod_{i^{\prime}\in \mathcal{N}(I_j)\setminus \{i,j\}}\Big(1-q\gamma_{i^{\prime}\rightarrow I_j}(1)\Big )\\
 & + \gamma_{j\rightarrow I_j}(0) \delta_{j\rightarrow I_j}(1)\Bigg[q-(1-q)\sum_{\ell=1}^{\vert \mathcal{N}(I_j) \vert -2}\sum_{S\in \Big\{\mathcal{N}(I_j)\setminus \{i,j\}\Big\}_\ell}\Big(-q\Big)^\ell\prod_{i^{\prime}\in S}\gamma_{i^{\prime}\rightarrow I_j}(1)\Bigg].
\end{split}\\
  \end{cases}
\end{equation}

In the end, when the algorithm is run for a fixed number $T$ of iterations, we compute the marginals of the posterior distribution for individuals at time $1$ as follows.
\begin{equation*}
\begin{cases}
{{\phi}}\Big(x_{i}^{(1)}=0\Big) \propto  \delta_{I_i\rightarrow i}(0) \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(0), \\
  {{\phi}}\Big(x_{i}^{(1)}=1\Big) \propto  \delta_{I_i\rightarrow i}(1) \displaystyle \prod_{t\in \mathcal{N}(i)}\mu_{t\rightarrow i}(1). \\
  \end{cases}
\end{equation*}
The LLRs of the marginals are given by
\begin{equation}\label{eq:llr1}
\lambda_i=\ln \frac{ {{\phi}}\Big(x_{i}^{(1)}=1\Big)}{{{\phi}}\Big(x_{i}^{(1)}=0\Big)}=\ln \frac{ \delta_{I_i\rightarrow i}(1)}{ \delta_{I_i\rightarrow i}(0)} + \sum_{t\in \mathcal{N}(i)}\ln \frac{\mu_{t\rightarrow i}(1)}{\mu_{t\rightarrow i}(0)}.
\end{equation}
The interpretation of the LLRs is done in the same way that has been explained in Section~\ref{subsec:bp}. For a given threshold $\tau$, individual $i$ at time $1$ is announced infected if $\lambda_i \geq \tau$.
Algorithm~\ref{alg:BP1} defines the belief propagation on combined graphs algorithm.
\begin{algorithm}
\caption{Belief Propagation on Combined Graphs}
\label{alg:BP1}
\label{CHalgorithm}
\begin{algorithmic}[1]
\State Initialize $\gamma_{i\rightarrow I_j}(1)=1-\gamma_{i\rightarrow I_j}(0)=p$ $\forall i\in [N], \forall j\in \mathcal{N}(i)$
\State Initialize $\mu_{t\rightarrow i}(0)=\mu_{t\rightarrow i}(1)=\frac{1}{2}$ $\forall t\in [M], \forall i\in \mathcal{N}(t)$
\For{$\ell=1,2,\cdots,T$ }
\State Compute $\delta_{I_i\rightarrow i}(0)$ and $\delta_{I_i\rightarrow i}(1)$ $\forall i\in [N]$ using~\eqref{eq:Intr-to-itm1}
\State Compute $\mu_{i\rightarrow t}(0)$ and $\mu_{i\rightarrow t}(1)$ $\forall i\in [N], \forall t\in \mathcal{N}(i)$ using~\eqref{eq:itm1-to-tst}
\State Compute $\mu_{t\rightarrow i}(0)$ and $\mu_{t\rightarrow i}(1)$ $\forall t\in [M], \forall i\in \mathcal{N}(t)$ using~\eqref{eq:tst-to-itm0} and~\eqref{eq:tst-to-itm1}
\State Compute $\delta_{i\rightarrow I_i}(0)$ and $\delta_{i\rightarrow I_i}(1)$ $\forall i\in [N]$ using~\eqref{eq:itm1-to-fctr}
\State Compute $\gamma_{I_j\rightarrow i}(0)$ and $\gamma_{I_j\rightarrow i}(1)$ $\forall j\in [N], \forall i\in \mathcal{N}(I_j)$ using~\eqref{eq:fctr-to-itm0j} and~\eqref{eq:fctr-to-itm0i}
\State Compute $\gamma_{i\rightarrow I_j}(0)$ and $\gamma_{i\rightarrow I_j}(1)$ $\forall i\in [N], \forall j\in \mathcal{N}(i)$ using~\eqref{eq:itm0-to-fctr}
\EndFor
\State Compute $\lambda_i$ $\forall i\in [N]$ using~\eqref{eq:llr1}
\end{algorithmic}
\end{algorithm}

\begin{example}\label{example}
\normalfont
Consider the combined graph shown in Fig.~\ref{fig:Comb}. We want to compute the BP messages exchanged over the edges of the combined graph. Calculating the messages from variable nodes to factor nodes is straightforward. Thus, we intend to compute messages from factor nodes to variable nodes. Since in~\cite{barbier2020strong,2018,coja2018belief,doi:10.1080/00018732.2016.1211393} it was shown that the messages from tests to individuals at time $1$ are computed using~\eqref{eq:tst-to-itm0} and~\eqref{eq:tst-to-itm1}, we only show how to compute messages from interaction nodes to individuals at time $0$ and time $1$.
The message from interaction node $I_j$ to individual $j$ at time $1$ is computed using
\[  \delta_{I_j\rightarrow j}\Big(\hat{{x}}_j^{(1)}\Big)   =\sum_{\hat{{x}}_i^{(0)}\in\{0,1\},~i\in \mathcal{N}(I_j)} \Bigg[P\Big({x}_j^{(1)}=\hat{{x}}_j^{(1)}\Big\vert \Big\{{x}_{i}^{(0)}=\hat{{x}}_{i}^{(0)}\Big\}_{i\in \mathcal{N}(I_j)}\Big)\prod_{i\in \mathcal{N}(I_j)}\gamma_{i\rightarrow I_j}\Big(\hat{{x}}_i^{(0)}\Big)\Bigg].
 \]
For instance, the message from interaction node $I_1$ to individual $1$ at time $1$ is given by
\begin{equation}\label{eq:intr1-to-itm1}
  \delta_{I_1\rightarrow 1}\Big(\hat{{x}}_1^{(1)}\Big)   =\sum_{\hat{{x}}_1^{(0)},\hat{{x}}_3^{(0)},\hat{{x}}_5^{(0)}\in\{0,1\}} P\Big({x}_1^{(1)}=\hat{{x}}_1^{(1)}\Big\vert {x}_{1}^{(0)}=\hat{{x}}_1^{(0)},{x}_{3}^{(0)}=\hat{{x}}_3^{(0)},{x}_{5}^{(0)}=\hat{{x}}_5^{(0)}\Big)\gamma_{1\rightarrow I_1}\Big(\hat{{x}}_1^{(0)}\Big)\gamma_{3\rightarrow I_1}\Big(\hat{{x}}_3^{(0)}\Big)\gamma_{5\rightarrow I_1}\Big(\hat{{x}}_5^{(0)}\Big).
\end{equation}
 We first consider the case that $\hat{{x}}_1^{(1)}=0$, and form Table~\ref{table:1}. It is easy to see that~\eqref{eq:intr1-to-itm1} can be expanded into the following
 \begin{table*}[h]
\centering
\begin{tabular}{|c|c|c|c|}
\hline
$\hat{{x}}_1^{(0)}$ & $\hat{{x}}_3^{(0)}$ & $\hat{{x}}_5^{(0)}$ &  $P\Big({x}_1^{(1)}=0\Big\vert {x}_{1}^{(0)}=\hat{{x}}_1^{(0)},{x}_{3}^{(0)}=\hat{{x}}_3^{(0)},{x}_{5}^{(0)}=\hat{{x}}_5^{(0)}\Big)$\\ \hline
0 & 0 & 0 & 1 \\ \hline
0 & 0 & 1 & $1-q$ \\ \hline
0 & 1 & 0 & $1-q$ \\ \hline
0 & 1 & 1 & $(1-q)^2$ \\ \hline
1 & 0 & 0 & 0 \\ \hline
1 & 0 & 1 & 0 \\ \hline
1 & 1 & 0 & 0 \\ \hline
1 & 1 & 1 & 0 \\ \hline
\end{tabular}
\caption{}\label{table:1}
\end{table*}

\[
\begin{split}
\delta_{I_1\rightarrow 1}(0) & =\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(0) +
(1-q)\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1)\\
& + (1-q)\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(0) + (1-q)^2\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1),
\end{split}
\]
where it can be simplified using the fact that messages $\gamma_{i\rightarrow I_j}$ are probability distributions on $\{0,1\}$, i.e., $\mu_{i\rightarrow I_j}(0) + \mu_{i\rightarrow I_j}(1)=1$.
\[
\delta_{I_1\rightarrow 1}(0)= \gamma_{1\rightarrow I_1}(0) \Big(1-q\gamma_{3\rightarrow I_1}(1)\Big)\Big(1-q\gamma_{5\rightarrow I_1}(1)\Big)
\]
We now consider the case that $\hat{{x}}_1^{(1)}=1$. Expansion of~\eqref{eq:intr1-to-itm1} results in
\iffalse
\begin{table*}[h]
\centering
\begin{tabular}{|c|c|c|c|}
\hline
$\hat{{x}}_1^{(0)}$ & $\hat{{x}}_3^{(0)}$ & $\hat{{x}}_5^{(0)}$ &  $P\Big({x}_1^{(1)}=1\vert {x}_{1}^{(0)}=\hat{{x}}_1^{(0)},{x}_{3}^{(0)}=\hat{{x}}_3^{(0)},{x}_{5}^{(0)}=\hat{{x}}_5^{(0)}\Big)$\\ \hline
0 & 0 & 0 & 0 \\ \hline
0 & 0 & 1 & $q$ \\ \hline
0 & 1 & 0 & $q$ \\ \hline
0 & 1 & 1 & $1-(1-q)^2$ \\ \hline
1 & 0 & 0 & 1 \\ \hline
1 & 0 & 1 & 1 \\ \hline
1 & 1 & 0 & 1 \\ \hline
1 & 1 & 1 & 1 \\ \hline
\end{tabular}
\end{table*}
According to this table, we have
\fi
\begin{multline*}
\delta_{I_1\rightarrow 1}(1)  =
q\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1) + q\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(0) + \Big(1-(1-q)^2\Big)\gamma_{1\rightarrow I_1}(0)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1)+\\
 \gamma_{1\rightarrow I_1}(1)\gamma_{3\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(0) + \gamma_{1\rightarrow I_1}(1)\gamma_{3\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1) + q\gamma_{1\rightarrow I_1}(1)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(0) + q\gamma_{1\rightarrow I_1}(1)\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1),
\end{multline*}
where we can simplify it to
\[
\delta_{I_1\rightarrow 1}(1)  = \gamma_{1\rightarrow I_1}(1) - \gamma_{1\rightarrow I_1}(0)\Big(-q\gamma_{3\rightarrow I_1}(1)-q\gamma_{5\rightarrow I_1}(1)+q^2\gamma_{3\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1) \Big).
\]

The messages from interaction nodes to individuals at time $0$ are computed as follows. The message from interaction node $I_j$ to individual $i$ at time $0$ is computed using
\[  \gamma_{I_j\rightarrow i}\Big(\hat{{x}}_i^{(0)}\Big)   =\sum_{\hat{{x}}_j^{(1)},\hat{{x}}_{i^{\prime}}^{(0)}\in\{0,1\}} P\Big({x}_j^{(1)}=\hat{{x}}_j^{(1)}\Big\vert {x}_{i}^{(0)}=\hat{{x}}_{i}^{(0)}, \Big\{{x}_{i^{\prime}}^{(0)}=\hat{{x}}_{i^{\prime}}^{(0)}\Big\}_{i^{\prime}\in \mathcal{N}(I_j)\setminus\{i\}}\Big)\delta_{j\rightarrow I_j}\Big(\hat{{x}}_{j}^{(1)}\Big)\prod_{i^{\prime}\in \mathcal{N}(I_j)\setminus\{i\}}\gamma_{i^{\prime}\rightarrow I_j}\Big(\hat{{x}}_{i^{\prime}}^{(0)}\Big).
 \]
For instance, the message from interaction node $I_1$ to individual $3$ at time $0$ is given by
\begin{equation}\label{eq:intr1-to-itm3}
 \gamma_{I_1\rightarrow 3}\Big(\hat{{x}}_{3}^{(0)}\Big)   =\sum_{\hat{{x}}_{1}^{(1)},\hat{{x}}_{1}^{(0)},\hat{{x}}_{5}^{(0)}\in\{0,1\}} P\Big({x}_1^{(1)}=\hat{{x}}_{1}^{(1)}\Big\vert {x}_{1}^{(0)}=\hat{{x}}_{1}^{(0)},{x}_{3}^{(0)}=\hat{{x}}_{3}^{(0)},{x}_{5}^{(0)}=\hat{{x}}_{5}^{(0)}\Big)\delta_{1\rightarrow I_1}\Big(\hat{{x}}_{1}^{(1)}\Big)\gamma_{1\rightarrow I_1}\Big(\hat{{x}}_{1}^{(0)}\Big)\gamma_{5\rightarrow I_1}\Big(\hat{{x}}_{5}^{(0)}\Big).
\end{equation}
First, we consider the case that $\hat{{x}}_{3}^{(0)}=0$, and expand~\eqref{eq:intr1-to-itm3} as follows.
\[
\begin{split}
\gamma_{I_1\rightarrow 3}(0) & = \delta_{1\rightarrow I_1}(0)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(0) + (1-q)\delta_{1\rightarrow I_1}(0)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1) + q\delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1)\\
& + \delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(0) + \delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1),
\end{split}
\]
where after simplification becomes
\[
\gamma_{I_1\rightarrow 3}(0)= \gamma_{1\rightarrow I_1}(1)\delta_{1\rightarrow I_1}(1)+ \gamma_{1\rightarrow I_1}(0)\delta_{1\rightarrow I_1}(0)\Big(1- q\gamma_{5\rightarrow I_1}(0)\Big)-\gamma_{1\rightarrow I_1}(0)\delta_{1\rightarrow I_1}(1)\Big(-q \gamma_{5\rightarrow I_1}(1)\Big).
\]
Then, we consider the case that $\hat{{x}}_{3}^{(0)}=1$.
\[
\begin{split}
\gamma_{I_1\rightarrow 3}(1) & = (1-q)\delta_{1\rightarrow I_1}(0)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(0) + (1-q)^2\delta_{1\rightarrow I_1}(0)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1) + q\delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(0)\\
& + \Big(1- (1-q)^2\Big)\delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(0)\gamma_{5\rightarrow I_1}(1) + \delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(0) + \delta_{1\rightarrow I_1}(1)\gamma_{1\rightarrow I_1}(1)\gamma_{5\rightarrow I_1}(1),
\end{split}
\]
where can be simplified to
\[
\begin{split}
\gamma_{I_1\rightarrow 3}(1) &= \gamma_{1\rightarrow I_1}(1)\delta_{1\rightarrow I_1}(1)+ (1-q)\gamma_{1\rightarrow I_1}(0)\delta_{1\rightarrow I_1}(0)\Big(1- q\gamma_{5\rightarrow I_1}(0)\Big)\\
&+\gamma_{1\rightarrow I_1}(0)\delta_{1\rightarrow I_1}(1)\Bigg[q-(1-q)\Big(-q \gamma_{5\rightarrow I_1}(1)\Big)\Bigg].
\end{split}
\]  \hfill \qedsymbol{}
\end{example}

\begin{figure*}[tb]
\begin{subfigure}{.49\textwidth}\centering \includegraphics[width=.95\linewidth]{sp1.eps}
  \caption{$\rho=0.01$}
  \label{fig:sub-second}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{sp5.eps}
  \caption{$\rho=0.05$}
  \label{fig:sub-fourth}
\end{subfigure}
\caption{Success probability as a function of the number of tests $M$ based on simulation results for $N=500$ individuals, the prevalence rate $p=0.01$, the contagion probability $q=0.1$, and the interaction probability $\theta=0.008$, under the binary symmetric noise model with parameter $\rho\in\{0.01,0.05\}$
.}
\label{fig:sp}\vspace{-0.25cm}
\end{figure*}

\section{Simulation Results}\label{sec:SR}
In this section, we compare the performance of the BP using Initial Prior probabilities (BPIP) algorithm, the BP using Updated Prior probabilities (BPUP) algorithm, and the BP on Combined Graphs (BPCG) algorithm using three metrics, success probability, FNR, and FPR.
Each result is averaged over $1000$ experiments. The testing matrix is constructed according to a Bernoulli design with parameters $\nu=\ln 2$. In the BPUP algorithm, the updated prior probabilities given by~\eqref{eq:prior} are computed using the contact tracing information. In our simulations, we assume that individual $i$ at time $0$, for each $i\in [N]$, interacts with individual $j$ at time $0$, for each $j\in [N]\setminus \{i\}$, with some fixed probability $\theta$, referred to as \emph{interaction probability}. It should be noted that $d_i$, the number of individuals which interact with individual $i$, follows a binomial distribution with parameters $N-1$ and $\theta$, i.e., $d_i\sim B(N-1,\theta)$. The expected value of $x_i^{(1)}$ is computed as follows:
\[\mathbb{E}\Big[x_i^{(1)}\Big]= \mathbb{E}\Big[ \mathbb{E}\Big[x_i^{(1)}\Big\vert d_i\Big] \Big] = \mathbb{E}\Big[1-(1-p)(1-pq)^{d_i}\Big]=1-(1-p)\mathbb{E}\Big[(1-pq)^{d_i}\Big],\]
where the term $\mathbb{E}\Big[(1-pq)^{d_i}\Big]$ is given by
\[
\mathbb{E}\Big[(1-pq)^{d_i}\Big]=\sum_{d=0}^{N-1}(1-pq)^d {N-1\choose d} \theta^{d}(1-\theta)^{N-d-1}=(1-pq\theta)^{N-1}.
\]
Thus, we have $\mathbb{E}\Big[x_i^{(1)}\Big]=1-(1-p)(1-pq\theta)^{N-1}$. Accordingly, the average number of infected individuals for the BPUP and the BPCG algorithms is given by $K=N\Big(1-(1-p)(1-pq\theta)^{N-1}\Big)$.

In Fig.~\ref{fig:sp}, we plot success probability as a function of the number of tests $M$ based on simulation results for $N=500$ individuals, the prevalence rate $p=0.01$, the contagion probability $q=0.1$, and the interaction probability $\theta=0.008$, under the binary symmetric noise model with parameter $\rho\in\{0.01,0.05\}$. The value of success probability for each number of test is optimized over the threshold in the range $\tau\in[-10,10]$. The number of iterations for the BPUP and the BPIP algorithms is $T=15$. We consider $T=30$ iterations for the BPCG algorithm.
It can be observed that the BPCG algorithm outperforms the other algorithms for all values of $M$. For instance, when $\rho=0.01$ and the number of test is $M=350$, the BPCG algorithm provides a success probability $4\%$ and $24\%$ greater than that of the BPUP and the BPIP algorithms, respectively.
Also, it can be seen that for the high noise regime, i.e., $\rho=0.05$, the advantage of BPCG algorithm over the other algorithms in terms of success probability becomes more evident. For example, for $\rho=0.05$ and $M=350$ tests, the success probability of the BPCG algorithm is $7.4\%$ and $43.5\%$ greater than that of the BPUP and the BPIP algorithms, respectively.

\begin{figure*}[tb] \vspace{-0.1cm}
\begin{subfigure}{.49\textwidth}\centering \includegraphics[width=.95\linewidth]{fnrfpr111.eps}
  \caption{$\rho=0.01$}
  \label{fig:sub-second}
\end{subfigure}
\begin{subfigure}{.49\textwidth} \vspace{-0.1cm}
  \centering
  \includegraphics[width=.95\linewidth]{fnrfpr55.eps}
  \caption{$\rho=0.05$}
  \label{fig:sub-fourth}
\end{subfigure}
\caption{FNR vs. FPR based on simulation results for threshold $\tau\in[-10,10]$, $N=500$ individuals, $M=350$ tests, the prevalence rate $p=0.01$, the contagion probability $q=0.1$, and the interaction probability $\theta=0.008$, under the binary symmetric noise model with parameter $\rho\in\{0.01,0.05\}$
.} \vspace{-0.25cm}
\label{fig:fnrfpr}
\end{figure*}

 In Fig.~\ref{fig:fnrfpr}, we depict FNR vs. FPR for all three decoding algorithms for threshold $\tau\in[-10,10]$, $N=500$ individuals, $M=350$ tests, the prevalence rate $p=0.01$, the contagion probability $q=0.1$, and the interaction probability $\theta=0.008$, under the binary symmetric noise model with parameter $\rho\in\{0.01,0.05\}$. Note that unlike the success probability, FNR and FPR do not converge. Instead, after a certain number of iterations, FNR and FPR oscillate around an average value. Thus, for each value of $\tau$, we compute the average FNR and FPR over a range of iterations. The range of iterations for the BPUP and the BPIP algorithms is $T\in\{15,16,\cdots,30\}$. For the BPCG algorithm, we consider the range of iterations $T\in\{30,31,\cdots,50\}$  .
 Each point on the curve corresponding to a decoding algorithm represents the pair (FNR,FPR) which has been computed for the same value of $\tau$. The closer a curve is to the origin of the FNR–FPR plane, the better the performance of the corresponding decoding algorithm in terms of FNR and FPR. It can be observed that for the BPCG algorithm the operating point that minimizes the total error rate, i.e., the sum of FPR and FNR, is closer to the origin than that of the BPUP and the BPIP algorithms.

\vspace{-0.17cm}

\bibliographystyle{IEEEtran}