A new COVID-19 detection method from human genome sequences using CpG island features and KNN classifier

Various viral epidemics have been detected such as the severe acute respiratory syndrome coronavirus and the Middle East respiratory syndrome coronavirus in the last two decades. The coronavirus disease 2019 (COVID-19) is a pandemic caused by a novel betacoronavirus called severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2). After the rapid spread of COVID-19, many researchers have investigated diagnosis and treatment for this terrifying disease quickly. Identifying COVID-19 from the other types of coronaviruses is a difficult problem due to their genetic similarity. In this study, we propose a new efficient COVID-19 detection method based on the K-nearest neighbors (KNN) classifier using the complete genome sequences of human coronaviruses in the dataset recorded in 2019 Novel Coronavirus Resource. We also describe two features based on CpG island that efficiently detect COVID-19 cases. Thus, genome sequences including approximately 30,000 nucleotides can be represented by only two real numbers. The KNN method is a simple and effective non-parametric technique for solving classification problems. However, performance of the KNN depends on the distance measure used. We perform 19 distance metrics investigated in five categories to improve the performance of the KNN algorithm. Some efficient performance parameters are computed to evaluate the proposed method. The proposed method achieves 98.4% precision, 99.2% recall, 98.8% F-measure, and 98.4% accuracy in a few seconds when any L1 type metric is used as a distance measure in the KNN.

SARS-CoV was first identified in the Guangdong province of southern China on 16 November 2002 [4] and MERS-CoV emerged in June 2012 in Jeddah, Saudi Arabia [5]. In late December 2019, a novel coronavirus SARS-CoV-2 emerged in China and spread rapidly throughout the world. The SARS-CoV-2 virus, which is genetically similar to SARS-CoV, caused a severe illness known as COVID-19 disease and a serious number of deaths worldwide. Since several early infected people visited a local seafood Market in December 2019 in Wuhan city of China, the virus is therefore thought to be a pathogen (SARS-CoV-2 pathogen) that jumped from an animal to a human and that caused an infectious disease. Similarity rate between SARS-CoV-2 and a bat coronavirus is 96% [6]. The outbreak COVID-19 was declared as a global pandemic on 11 March 2020 by the World Health Organization (WHO). As of 8 October 2020, the disease spread to 188 countries and territories has infected over 36.2 million people and has caused more than 1.05 million deaths; more than 25.2 million people have recovered from this illness.
The most common symptoms of COVID-19 disease are cough, fever, gastrointestinal and musculoskeletal symptoms as well as loss of taste or smell. Shortness of breath and chest pressure or pain are among less common symptoms [7]. Because these common symptoms are similar to the common flu, it is difficult to conduct an early diagnosis of SARS-CoV-2. It is essential to quickly identify positive cases since it spreads rapidly and poses a threat to the public health system. Furthermore, there is still no specific antiviral drug recommended for the treatment of the novel coronavirus disease other than supportive care so far. In the period of pre-clinical researches, some medications already prescribed for other diseases have shown positive effects against COVID-19 virus. Due to the absence of an effective COVID-19 treatment, early detection of this life-threatening infectious disease plays a very important role in all medical therapies and in the prevention of the spread of the disease. Recently, many data scientists have been working extensively on remarkable features of the virus. For this purpose, artificial intelligence applications and machine learning methods have also been used successfully to speed up the diagnosis process of the COVID-19 cases [8][9][10]. Therefore, the learning based classification models will not only reduce the burden on healthcare professionals but will also facilitate early diagnosis of the disease.
In this paper, we introduce an effective classification method to distinguish SARS-CoV-2 from common human coronaviruses, which are AlphaCoV, BetaCoV-1, MERS-CoV, HKU1-CoV, NL63-CoV, and 229E-CoV using CpG island features and KNN classification method. Main contributions of this research can be summarized as follows: The choice of the differentiable features based on characteristic of SARS-CoV-2 is a critical step to improve classification performance. In this study, we propose CpG island features to differentiate SARS-CoV-2 from the other human coronaviruses. We propose a robust prediction method by using KNN classifier with any L 1 type distance metric selected from among 19 distance metrics investigated in five categories, which are L 1 type, L 2 type, vicissitude, inner product metrics, and the other types of metrics. We construct a larger dataset containing almost all types of these human coronaviruses which are genetically similar to SARS-CoV-2. They are AlphaCoV, BetaCoV_1, MERS-CoV, NL63-CoV, HKU1-CoV, and 229E-CoV. The proposed method efficiently detects COVID-19 cases in a few seconds on the relatively large dataset.
The rest of the paper is organized as follows. In Section 2, a literature survey on COVID-19 is presented. In Section 3, the Knearest neighbors (KNN) method and distance measures used in the KNN method are summarized. The proposed classification strategy to detect COVID-19 cases is introduced in Section 4. Section 5 reports experimental results and evaluation of the experimental observations conducted by various distance metrics. In Section 6, we compare the proposed method with the other detection methods. Finally, Section 7 presents concluding remarks and future directions.

Related works
Deep learning and machine learning techniques have been widely used in a variety of research areas such as big data analysis [11][12][13], image classification [14], face detection [15,16], and disease prediction [17][18][19][20][21]. The computer based diagnostic systems developed by the help of artificial intelligence techniques will speed up early diagnosis of COVID-19 and thus will decrease the workload on health-workers. Chen et al. [22] published a comprehensive survey about studies performed by using artificial intelligence in the literature on the fight against COVID-19. In addition, they observed that machine learning, deep learning and artificial neural networks technologies have been successfully implemented at almost every stage of combating COVID-19 when compared to other coronaviruses such as SARS-CoV and MERS-CoV.
Randhawa et al. [23] suggested a combination of supervised machine learning with digital signal processing methods (ML-DSP) for an accurate and scalable taxonomic classification of genomic sequences. They matched each genomic sequence to discrete values corresponding to its genomic signals. They employed six supervised machine learning classifiers (linear discriminant, linear support vector machine, quadratic support vector machine, fine KNN, subspace discriminant, and subspace KNN) to detect SARS-CoV-2. They used 29 SARS-CoV-2 genome sequences and 20 genome sequences for each of alphacoronavirus, betacoronavirus, and delta-coronavirus. They observed that the linear discriminant method achieved 100% accuracy. Naeem et al. [24] proposed a method to distinguish COVID-19, SARS-CoV, and MERS-CoV viruses by using the K-nearest neighbors and the trainable cascadeforward back-propagation neural network methods. They extracted genomic signal processing features using a dataset that contains 76 genome sequences for each type of coronavirus from the National Center for Biotechnology Information (NCBI). Their results showed that performance of the KNN algorithm was higher than the cascade-forward back-propagation neural network in all COVID-19/SARS-CoV, COVID-19/MERS-CoV and COVID-19/SARS-CoV/MERS-CoV classification processes, and achieved an accuracy of 100%. Batista et al. [25] developed a new method to predict whether patients in the emergency care unit are COVID-19 positive or not by using five supervised machine learning algorithms, which are neural networks, random forests, logistic regression, support vector machines, and gradient boosting regression trees. They collected data from 235 patients, and they used 15 different types of features such as age, gender, hemoglobin etc. Their experimental results showed that the support vector machine classifier achieved the highest performance with AUC value of 84.7% among the other classifiers.
Unal and Dudak [26] studied on diagnosis of COVID-19 viral disease. They applied Naive Bayes, KNN, support vector machines and decision tree algorithms to the dataset named as COVID-19 Mexico Patient Health Dataset. The dataset consists of 95839 cases recorded by Mexican government, and 19 different types of features like as the sex of the patient, age of the patient, the state of pneumonia and intubation as well as the state of many other diseases. They performed four types of supervised machine learning algorithms. Their experimental results showed that support vector machine achieved the best predictive performance with the classification accuracy of 100%.
Although there are only few studies detecting COVID-19 cases from genome sequences, a number of papers have been published for detecting COVID-19 cases from X-ray or computed tomography (CT) images [27][28][29][30][31][32][33][34], recently. Barstugan et al. [27] developed a COVID-19 classification method by using the support vector machine classifier on a dataset including different types of 150 abdominal CT images. Their experimental results showed that their proposed method achieved an accuracy of 99.6%. Ozturk et al. [28] applied the DarkCovidNet deep learning model to raw chest X-ray images and generated a binary classification (COVID-19, no findings). They also used the DarkCovidNet classifier method on the same dataset to create a triple classification (COVID-19, nofindings, pneumonia). They declared that the highest accuracy of the classifier was 98.08% for the binary classification and 87.02% for the triple case. Sekeroglu et al. [29] proposed an alternative COVID-19 detection method by using deep learning and machine learning classifiers on a publicly available dataset which contains 1583 healthy, 4292 pneumonia and 225 confirmed COVID-19 chest X-ray images. They stated that a convolutional neural network (CNN) without pre-processing and with minimized layers achieved an accuracy of 98.50%. Jain et al. [30] developed a new approach to detect COVID-19 cases among chest X-ray images of healthy, bacterial pneumonia, viral pneumonia and COVID-19 by using ResNet18, ResNet101, DenceNet121 and VGG-16 deep learning models. Their experimental results demonstrated that ResNet101 method achieved the highest accuracy of 98.93%. Apostolopoulos and Bessiana [31] used VGG19, Mobile Net, Inception, Xception and Inception ResNet v2 deep transfer learning classifiers to detect COVID-19 positive cases among 224 COVID-19, 700 bacterial pneumonia and 504 normal chest X-ray images. Their results presented that VGG19 achieved the best binary classification accuracy of 98.75% over the other CNN methods.
Ahuja et al. [32] used deep learning methods to detect COVID-19 positive cases from chest X-ray images. Asnaoui and Chawki [33] constructed a novel method for detecting COVID-19 from pneumonia chest X-ray and tomography images. They performed recent deep learning methods, VGG16, VGG19, DenseNet201, Inception_ResNet_V2, Inception_V3, Resnet50, and MobileNet_V2. They reported that the performance of Inception_Resnet_V2 demonstrated the best accuracy of 92.18%. Basu et al. [34] proposed an alternative screening method of COVID-19, which is called Domain Extension Transfer Learning. They extracted some discriminate features from the chest X-ray dataset. They then classified the images as normal, pneumonia, other diseases, and COVID-19 with an accuracy of 95.3%.

Background
We propose a new COVID-19 detection method based on CpG island features and the KNN classifier. The KNN is extremely useful for large data classification and its performance depends on the distance metric used [35][36][37][38]. Several studies have been conducted to detect optimum metrics for KNN algorithm in [39,40]. In this section, first, we provide a brief description of the KNN method and some recent studies based on the KNN. Second, we investigate the metrics used in the KNN classifier under five categories to improve the performance of the KNN classifier on our model.

K-nearest neighbors
K-nearest neighbors (KNN), a supervised machine learning algorithm, can be efficiently used to solve classification problems. The KNN was introduced in 1951 by [41] and then recasted in 1967 by [42]. The KNN is a non-parametric classifier known as one of the simplest and laziest algorithms. That is, there is no need to create a learning model in this classification method. Despite the lazy structure of KNN, it was proposed as one of the 10 most effective methods in the process of analyzing information in a database given [43]. In the prediction process, the class to which a new observation data belongs is determined by calculating the shortest distance between the observation sample and its K-nearest neighbors samples.
There are some recent studies to improve the performance of the KNN. To decrease the sensitivity of the neighborhood size of k and improve voting strategy in the region of neighborhoods, Gou et al. [44] proposed two k-nearest neighbor rules, which are the weighted representation-based k-nearest neighbor rule and the weighted local mean representation-based k-nearest neighbor rule. Their experimental results demonstrate that the proposed methods have a lower sensitiveness to k. Gou et al. [45] proposed another study to improve the selection of the neighborhood size of k by introducing the generalized mean distance-based KNN classifier. They stated that the proposed method is less sensitive to k over the KNN-based classifier.

Metrics
To improve the performance of the KNN classifier, we investigate the metrics used in the KNN algorithm under five categories by following the similar categorization in [39,40]. We list them as below:

Proposed method
In this section, we present a new COVID-19 detection method based on the KNN classifier and CpG island features. Main steps of the proposed method are described in 1. The first step of the algorithm is feature extraction. Extracting robust and discriminative features from human coronaviruses genome sequences is the most critical step to improve diagnosis of SARS-CoV-2. In this step, we propose to use CpG based features since CpG dinucleotides in the open reading frames of SARS-CoV-2 has extremely lowfrequency [46]. The main reason of lower CpG dinucleotides density is the mutation of C into A and G into T [47,48]. The proposed features are extracted by using Eq. 1 and Eq. 2.
CGp ¼ ratioðCÞ þ ratioðGÞ ð 1Þ where ratioðCÞ; ratioðGÞ, and ratioðCGÞ are computed as divided the number of occurrences of C; G, and CG in the sequence, respectively by the sequence length. Thus, each sequence containing 30 000 nucleotides is represented by two features only. 1 provides an example of how features are calculated from a part of the sequence. After feature extraction step, we apply the KNN method to classify SARS-CoV-2 sequences. The performance of the KNN mainly depends on the metric that is used to compute the distances between different data samples. To improve the performance of the KNN algorithm, we perform 19 distance metrics investigated in five categories, which are L 1 type, L 2 type, vicissitude, inner product metrics, and the other types of metrics. In this step, we propose to use L 1 type metric as a distance measure in the KNN.

Algorithm 1. Proposed Algorithm
Require: S: training genome sequences, L: COVID-19 positive, COVID-19 negative seq: a test sequence, and k: the neighborhood size Ensure: Determine the class label of the test sequence seq Step 1: Feature Extraction 1: for all genome sequences do 2: feature 1 ¼ ratioðCÞ þ ratioðGÞ 3: feature 2 ¼ ratioðCGÞ=ratioðCÞratioðGÞ 4: end for Step 2: Apply KNN and use L 1 type metric 5: Compute the distance between seq and every sample in S using any L 1 type metric 6: Choose k samples in S that are nearest to seq 7: Assign seq to majority class

Experimental results
In this section, we evaluate performance of the proposed method. Before discussing the results, first we explain our dataset, second we mention experimental setup, and then we provide a brief information about the performance measures. Finally, we will give the results of the experimental observations conducted by using different distance metrics.

Dataset
The 2019 novel coronavirus resource (2019nCoVR) [49] at China National Center for Bioinformation is one of the most important sources of various types of coronaviruses. It integrates various important databases including the GISAID, NCBI, NMDC and CNCB/NGDC. In this study, we used complete genomic sequences of human coronaviruses obtained from 2019nCoVR. Genome sequence of each coronavirus has approximately a length of 30 000 nucleotides. The properties of human coronavirus sequences are presented in 6. The dataset includes various types of human coronaviruses such as AlphaCoV, BetaCov-1, MERS-CoV, NL63-CoV, HKU1-CoV and 229E-CoV as well as SARS-COV-2. We refer SARS-CoV-2 sequences as COVID-19 positive, and the sequences that do not include SARS-CoV-2 are referred as COVID-19 negative. In addition to 1000 SARS-CoV-2 sequences, we used 592 genome sequences of other human coronaviruses in our experiments. We note that all available genome sequences of human coronaviruses, which are different from SARS-CoV-2 are downloaded.

Experimental setup
The experiments were performed using a core i7, 2.7 GHz processor, 16 GB RAM under Linux operating system. Feature extrac-

Table 6
The properties of complete genome sequences of human coronaviruses. tion and the classification process were performed using Python programming language. Precision, recall, F-measure, and accuracy in predicting COVID-19 positive case are used as the major performance measures. Next, we define these performance measures.

Performance measures
COVID-19 prediction as a binary classification problem has four prediction outcomes. True Positive (TP) is the number of genomic sequences which are correctly classified as COVID-19 positive, True Negative (TN) is the number of genome sequences which are correctly classified as COVID-19 negative, False Positive (FP) is the number of genome sequences which are incorrectly classified as COVID-19 positive, and False Negative (FN) is the number of sequences which are incorrectly classified as COVID-19 negative. The performance evaluation metrics such as precision, recall, Fmeasure, and accuracy are defined by using these outcomes and presented in confusion matrix by assigning actual and predicted labels in 2. Next, we briefly explain performance parameters used in this study.
Precision is the accuracy of the classifier in the presence of false positive case. It is computed as the ratio of the number of correctly classified positive samples to the number of samples labeled by the system as positive and it is exhibited in Eq. 3.
Recall refers to as the number of positive class predictions made out all positive samples in the data set. It is computed as in Eq. 4.
F-Measure is determined by the harmonic mean of precision and recall, and computed using Eq. 5.
Accuracy is computed by dividing the total number of true cases by all cases and it is indicated in Eq. 6.

Evaluation of the proposed features and KNN classifier with respect to metrics
In this section, we evaluate the efficiency of proposed features and the performance of the KNN classifier with respect to metrics using a fivefold cross-validation that was performed on the dataset. The dataset is randomly divided into two sets which are training and testing. Eighty percent of the entire human genome sequences is used for training and the remaining 20% of the dataset is used for validation. The experiments are repeated five times as presented in 3. In our experiments, we increase the neighborhood size of k between 1 and 20, and observe that the classification performance of the proposed method remains the same or slightly decreases with increasing the neighborhood size of k. For this reason, in all experiments, the k value is basically set to 1.
7 presents precision, recall, and F-measure values of the KNN classifier with respect to five metric groups, which are L 1 type, L 2 type, vicissitude, inner product metrics, and the other types of metrics. The results of the metrics in the same group are close to each other; thus, we present the average of the results of the experiments performed by the metrics in the same group. First, we analyze the results of the KNN equipped with L1 type metrics, which are KM, ChebM, ManM, SM, CanM, and MCM. The KNN with L1 type metrics presents the best result, and achieves a precision of 98.4%, a recall of 99.2%, and an F-measure of 98.8%. This means that almost all sequences are classified correctly. Next we look at the results of the KNN with L2 type metrics, which are ClaM, EM, SquM, NCSM, DivM, and SCSM. The KNN classifier with this group metrics achieves a precision of 96.0%, a recall of 98.2%, and an Fmeasure of 97.1% on average. When we look at the results of the KNN with vicissitude metrics, it presents better result than KNN with L2 type metrics. It achieves a precision of 98.4%, a recall of 99.0%, and an F-measure of 98.7%. Next, we investigate the results of the KNN with inner product metrics, ChoM and DicM. This group of metrics presents the worst results among the other groups, and achieves a precision of 94.4%, a recall of 92.8%, and an F-measure of 93.4%. Finally, we investigate the other type of metrics, which are HasM and MotM. The results of KNN with HasM are close to the results of the KNN with inner product type metrics, and achieves a precision of 94.4%, a recall of 92.9%, and an F-measure of 93.4%. On the other hand, the results of the KNN with MotM is close to them of L2 type metrics, and are better than that of HasM. It achieves a precision of 96.4%, a recall of 98.2%, and an F-measure of 97.3%.
In addition to precision, recall, and F-measure values, we present accuracy results. We take the average of the accuracies of KNN classifier obtained with the metrics in the same group, and show the results for each metric group separately in 4. The accuracy results of the KNN classifier using inner product metrics and HasM are close to each other and have the worst accuracy, 91.7%, and 91.8%, respectively. When L2 type metrics are used, the method achieves an accuracy of 96.2%. The KNN with MotM metric has remarkable accuracy values, which is 96.5%. The accuracy values of the KNN with L1 and Vicissitude type metrics are close to each other, and our method achieves the best accuracy with 98.4% when L1 type metrics are used.

Comparisons with other methods
In this part, we compare the results of the proposed method with the state-of-the-art methods. 8 provides a simple comparison of the results in terms of the method, dataset, class, and accuracy. We investigate these studies under two categories: the studies using genome sequences datasets and the other datasets.
First, we compare our method with the studies which detect COVID-19 cases from genome sequences. Randhawa et al. [23] combined supervised machine learning methods with digital signal processing. They used 29 sequences of SARS-CoV-2, and 20 sequences for each of alphacoronavirus, betacoronavirus, and deltacoronavirus. They achieved 100% accuracy when linear discriminant method was used. Main limitation of their proposed method is the number of sequences used in their dataset. Furthermore, Randhawa et al. used delta coronaviruses genomes. However, deltacoronavirus mainly causes an infectious disease among bird species rather than human. When the number of sequences in their dataset is increased and human coronavirus sequences genetically similar to SARS-CoV-2 sequences are added to their dataset, overall accuracy of their method may decrease. Another study predicting COVID-19 from genome sequences is introduced by Naeem et al. [24]. They used the classical KNN method to distinguish SARS-CoV-2 sequences among SARS-CoV-2, SARS-CoV and MERS-CoV genome sequences. They used 76 sequences for each of SARS-CoV-2, SARS-CoV, and MERS-CoV, and achieved an accuracy of 100%. They extracted features by using Discrete Fourier transform, Discrete cosine transform, and Seven Moment Invariants methods. When comparing to our method, their feature    [24] is that they work with a smaller dataset just like Randhawa.
Next we discuss advantages/drawbacks of the methods using the other datasets. Unal and Dudak [26] used Mexico Patient Health Dataset to detect COVID-19 cases. Although the size of their dataset is large, their dataset only represents the features of a specific region. Furthermore, they used 19 different features including the sex of the patient, age of the patient, the state of pneumonia and intubation as well as the state of many other diseases. One of the advantages of our model over the method of Unal and Dudak [26] is to use powerful and effective two features derived from the complete genome sequences of human coronaviruses. Thus, our proposed method detects COVID-19 positive cases within a few seconds. However, Unal and Dudak [26] did not provide any explanation about the detection time. Another different aspects of the proposed method from their method is that they use the classicial KNN classifier with default parameters although we use the KNN classifier with an optimum distance metric in our proposed model.
There are a number of the existing studies detecting COVID-19 cases from image dataset in the literature. In 8, we also exhibit six image based studies achieving remarkable accuracy values. These image based methods used different deep learning or machine learning methods. Although image based studies have a higher accuracy, when considering a high mutation rate of SARS-CoV-2, using genomic sequences is extremely beneficial when tracking coronavirus genes that change frequently as the disease spreads from one person to another. Moreover, the radiation generated by X-ray or CT scanning machines may cause permanent damages to people. For this reason, X-ray or CT scanning may not be obtained for some people, which can be considered as a disadvantage of image based studies.

Conclusion
COVID-19 is an existing epidemic that sets new records in terms of cumulative and daily numbers for global infection. The pandemic is an unprecedented situation in healthcare systems worldwide, and to overcome this pandemic, it is essential to accurately detect COVID-19 cases by analyzing data of patients in a minimum amount of time. In this study, we propose an accurate and fast method to detect COVID-19 positive cases from genome sequences of human coronaviruses. In the proposed method, first, the features that significantly differentiate COVID-19 cases are extracted from the complete genome sequences of human coronaviruses. In this step, we propose to use CpG island features. Each genome sequence of human coronavirus, which includes about 30,000 nucleotides is represented by two real numbers only. Feature extraction step takes just a few seconds. Second, the KNN method is used for the classification of COVID-19 positive cases from the other types of human coronaviruses, AlphaCoV, Beta-Cov-1, MERS-CoV, NL63-CoV, HKU1-CoV, and 229E-CoV. The KNN classifier is the simplest method and has high flexibility for solving complex classification problems. The accuracy of the KNN is higher than state-of the-art classifiers in certain cases, and it often produces efficient performance. However, the performance of the KNN greatly depends on the metric performed. To detect the most appropriate metric, we review five groups of metrics used in the KNN classifier. The selection of different distance metrics for the KNN can result a variation in accuracy outcomes for the same dataset. Experimental results reveal that the proposed method achieves the highest accuracy, which is 98.4% on average in a few seconds when L 1 type metrics are used as a distance measure in the KNN. In future studies, we will compare human SARS-CoV-2 sequences to other types of coronavirus sequences such as bat SARS-CoVlike coronaviruses 2, and propose a similarity based feature to increase overall accuracy. In addition, the factors affecting the recovery status of patients suffering from COVID-19 may be investigated in future studies by combining machine learning and parallel computing methods with effective features.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


qprof: a gprof-inspired quantum profiler

We introduce qprof, a new and extensible quantum program profiler able to generate profiling reports of various quantum circuits. We describe the internal structure and working of qprof and provide three practical examples on practical quantum circuits with increasing complexity. This tool will allow researchers to visualise their quantum implementation in a different way and reliably localise the bottlenecks for efficient code optimisation.

I. INTRODUCTION
The quantum computing field has been evolving at an increasing rate in the past few years and is currently gaining more traction. Several quantum chips, the underlying hardware that enable researchers and companies to run quantum algorithms, have been announced by different research teams. The error rates and number of qubits provided by these chips greatly improved, with quantum hardware that have up to 65 qubits in early 2021.
Software has also seen a tremendous rise with the emergence of several quantum computing frameworks and languages such as Qiskit [2], Q# [29], PyQuil [7], Cirq [6] or myQLM [28] to name a few. These frameworks help in speeding-up the process of implementing a quantum algorithm by providing their own "standard library". Most of them also include specialised libraries whose purpose is to facilitate the development and testing of new quantum algorithms. For example, all the quantum computing frameworks cited previously include a library to simulate quantum circuits, some even implement several simulation algorithms such as a full state-vector simulator, a simulator for stabilizer circuits [1,15] or a simulator using matrix-product states as described in [25,30]. Most of the frameworks that target real quantum chips also include libraries to characterise a given quantum hardware, using for example randomised benchmarking [8,10,14,20,23] methods, or to mitigate hardware noise [4,21].
Finally, a large majority of the quantum computing frameworks provide a way to automatically optimise a quantum circuit. This optimisation is often performed during compilation, when the abstract quantum circuit representation is translated to be compliant with the targeted hardware. Automatic optimisation of quantum circuits is a broad area of research with algorithms based on pattern-matching [19,22,24], gate-based optimisation algorithms [3,13] or even pulse-based ones [9,16,26].
But even though automatic optimisation has already been shown to be successful in optimising complex quantum circuits [5], most algorithms only perform local op- * adrien.suau@lirmm.fr timisations without any knowledge about the algorithms used and do not have a global vision of the optimised circuit. This means that automatic optimisation algorithms will likely not be able to perform some optimisations that require a global vision of the quantum circuit or a semantic understanding of the quantum algorithms used. Identifying the usage of a non-optimal algorithm in the implementation and replacing it with a more efficient algorithm is an example of such an optimisation that cannot be performed without a semantic understanding of the implementation.
Currently, the only other way one has to optimise a given quantum implementation is "trial and error": try to locate a "hot spot" in the implementation, either by a tedious theoretical analysis or a manual counting of the routine calls, optimise the potential hot spot found and finally check if the optimisation performed improved the performances. This process has a severe drawback that makes it impractical on real-world implementations: the first step that consists in finding the hot spots is either imprecise or potentially very long, tedious and errorprone on large implementations. qprof aims at replacing this manual, tedious and errorprone step of analysis by automatically generating a report with all the useful information needed to find the hot spots of the given quantum program implementation. The qprof tool has been strongly inspired by classical profilers such as gprof [11,17] which try to solve the exact same issue but in classical (non-quantum) programming.
The rest of the paper is organised as follows. In Section II we explain the internals of qprof and explain in details its architecture, the design choices made, and their impact on the tool efficiency, extensibility and usability. We then include code snippets and simple practical examples in Section III to illustrate the tool usage.

II. HOW DOES QPROF WORKS?
A. General structure The general structure of qprof is composed of 3 main parts that interact with each other: framework plugins, core data structures and logic, and exporters.

arXiv:2106.07326v1 [quant-ph] 14 Jun 2021
The overall workflow of qprof is schematically explained in Figure 1. In this workflow, qprof can be seen as a black-box that takes a "quantum circuit" as input and returns a "profiler report". This black-box view should be enough for users that only want to use the qprof tool, but experienced users or plugins developers might need more details on the internals of qprof in order to understand how it works.
The following sections will introduce in details the three different parts that compose qprof. Section II B describes the plugins architecture, used by qprof to be as modular and extensible as possible. Section II C then shows how this plugin architecture is used to allow qprof to be mostly framework-agnostic by allowing anyone to implement a framework without modifying qprof code. A description of the core data structures and core logic is then provided in Section II D. Finally an explanation of the different exporters natively provided by qprof is given in Section II E.

B. Plugin architecture
The qprof tool aims at being the standard for profiling quantum circuits, independently of the framework they are written with. In order to be versatile and support as much current and future quantum computing frameworks as possible, qprof has been designed to be easily extensible to support new frameworks. qprof extensibility is achieved via a system of runtimediscovered plugins. In order to be discoverable by qprof, a plugin should naturally obey some rules such as implementing a specific interface (explained in Section II C) or being organised in an imposed way. But the key ingredient to extensiblity is that plugins do not have to be part of the main qprof tool: they can be developed and published by anyone. This allows several situations that may help improving qprof stability and evolution in time. For example, users might decide to roll-out their own plugin to support a new framework they are using internally. Another important situation that is made possible by this plugin architecture is that framework vendors have the opportunity to provide a qprof plugin along with their framework and to maintain it as an official plugin.
Finally, such an architecture based on external plugins allows the user to only install the plugins and frameworks needed instead of installing all the plugins along with qprof. This simple side improvement greatly reduces the installation time, installation size and plugin discovery time as it avoids installing and loading unused quantum computing frameworks. Observation 1. Even though the reduction in installation time, size and loading time may seem negligible at first sight, we found out that avoiding some packages adds a significant improvement. For example, we computed that installing qiskit 0.24.1 with Python 3.8.8 on a modest laptop (MacBook Air from 2017) with a very good Internet connection takes more than 1 minute and ends up using more than 600 Mo of disk space. About runtime effects, importing the main qiskit package from a Python interpreter can take up to several seconds.

C. Framework support
The goal of the qprof tool is to provide a unique profiling interface for quantum programs. As such, it should be able to support as much quantum computing programming frameworks as possible. Taking into account that several of the most successful quantum computing frameworks such as Qiskit, Cirq, PyQuil or myQLM are Python libraries, and in order to ease its integration with these already existing frameworks, qprof has naturally been designed as a Python library too.
But being written in the same language as a given quantum computing framework does not grant a direct compatibility. In fact, all the previously cited frameworks use a different representation for quantum programs. Qiskit implements a QuantumCircuit class, Cirq comes with a Circuit class, PyQuil uses a Program structure and finally myQLM uses either QRoutine or AbstractGate instances.
To be as independent as possible from the different quantum computing frameworks and their internal representation of a quantum circuit, qprof has been designed to be both extensible and generic.
Extensibility is achieved with the help of the plugin architecture described in Section II B. And in order to be as generic as possible, qprof uses an abstract common interface to represent the concept of "quantum (sub)routine". This concept is formally defined in Definitions 1 to 3.

Definition 1. Quantum routine:
a (possibly parametrised) named sequence of quantum subroutines. Definition 2. Quantum subroutine: a quantum routine that is part of a higher-level quantum routine (i.e. that is called by another quantum routine).
Definition 3. Native quantum subroutine: a quantum subroutine that represents a native hardware operation and that does not call any quantum subroutine.
Using Definitions 1 to 3, a common interface for the concept of "quantum routine" emerges. First, a quantum routine should have a name that can be retrieved. Secondly, we should be able to distinguish between native quantum subroutines and non-native ones. Finally, for non-native quantum subroutines, we need a way to iterate over all the subroutines composing it.
This interface, schematised in Figure 2, is the core abstraction layer of qprof that allows it to be as independent as possible from the underlying quantum computing framework used to represent the profiled quantum circuit. Each framework-specific representation is wrapped by a RoutineWrapper. Each supported framework should have a corresponding qprof plugin that implements the RoutineWrapper interface. The init method initialises a RoutineWrapper instance with an instance of the framework-specific quantum circuit representation. The name method returns the name of the currently wrapped routine. is base returns True if the routine is a native subroutine as defined in Definition 3, else it returns False. Finally, the iter method returns an object that can be iterated on and whose iterates are the different subroutines called by the current subroutine.

D. Core data structures and logic
Once the issue of adapting qprof to the various quantum computing frameworks has been solved, we can start considering the main problem of profiling a quantum circuit. The first important question to answer is: what kind of figures are interesting to include in a profile report? Section II D 1 tries to address this question as extensively as possible. Then, once the target quantity to measure is known, our main concern is: how to compute this quantity from the abstract quantum circuit representation presented in Section II C? qprof approach to solve this issue is presented in Section II D 2.

Interesting data to profile
Profiling a program is the action of gathering data on its execution. For classical programs and profilers, the list of data that can be gathered is quite extensive ranging from high-level quantities such as the time spent in a given function or the memory used during the program execution to low-level information recovered via hardware counters such as cache misses or branch-predictionmisses.
But for quantum computing, the quantities of interest need to be adapted as several classical data such as cachemiss or branch-prediction-miss do not have any meaning anymore. Nevertheless some classical quantities have a quantum analogue that may be useful for optimisation purposes. This is the case for the classical "instruction number" quantity, that translates trivially to its quantum counterpart "native gate number" (or "hardware gate number"). The number of native gates executed by a quantum routine is a useful information for several reasons: it is simple, the routine worst-case execution time can be computed from it and a lower-bound of the routine error rate can also be devised using this information.
Another classical quantity that can be translated to quantum computing is the "time spent in routine". This quantity can be subdivided in two more specific figures: the "time spent exclusively in routine" (sometimes called "self time") and the "time spent in subroutines called by the routine". This separation is often done in classical profiling programs as having these two execution times gives very useful information about the profiled routine that cannot be obtained from the "time spent in routine" only.
The last classical quantity with a meaningful quantum counterpart is the "memory usage", which may be translated as "number of qubits needed" when using quantum computers.
About quantities without a clear classical parallel but potentially useful, one can cite the "routine depth" as an approximation of the total execution time of the routine, the "T-count" for error-correction estimates or the "idle time" to estimate the potential effects of qubit decoherence on the routine.

Graph representation (call-graph)
Following Definitions 1 to 3 and the RoutineWrapper interface we defined in Section II C, a graph-like representation of a quantum program seems to be particularly well suited. In this representation, nodes are quantum subroutines and an oriented edge from node A to node B means that the quantum subroutine represented by A calls the quantum subroutine represented by B. This representation of a program is called a call-graph in classical computing.
Using Definitions 1 to 3 and the RoutineWrapper in- init , name, is base and iter methods are described in Figure 2. The eq method tests if other routine is equal to the current instance. hash computes an integer hash of the currently wrapper routine.
Grover's algorithm Call-graph representation of one possible implementation of Grover's algorithm. Each node represents one routine rather than one call to the routine. Edges have been regrouped and labelled with the number of calls for readability purposes. In the internal representation used by qprof, edges are not regrouped and are ordered to account for the original quantum program subroutine call order. terface, Figure 3 shows a representation of one possible Grover's algorithm implementation. Even though this representation is valid according to the general definition of a call-graph, it contains a lot of redundant information that scrambles the useful data in visual noise. Because of this, most of the call graphs representations avoid the duplication of nodes, i.e. create one node for a specific subroutine and re-use this node whenever the subroutine is called. And this brings to the main optimisation performed by qprof under-the-hood to keep its efficiency even on very large quantum circuits: routine caching. Instead of blindly performing a full exploration of the call-graph, qprof will keep in memory all the previously encountered subroutines and re-use the already computed data whenever possible. In order to make the cache efficient, qprof uses a hash-map structure. This means that in order to use the cache mechanism, RoutineWrapper instances should now be also hashable and comparable for equality. The final RoutineWrapper interface that includes these requirements is represented in Figure 4.
Using such an optimisation, the call-graph represented in Figure 3 changes to the one in Figure 5.
It is important to realise that this optimisation does BaseExporter init (**args) export(ProgramData) -> Any FIG. 6. Exporter interface. Any plugin that implements an exporter should be a derived from the BaseExporter class and implement this interface. The return type of the export method is not specified and can be anything. The main profile function will return the output of the exporter.
not delete any kind of information when compared to the non-optimised version. As noted in Figure 5, the order of subroutine calls is preserved and this optimisation simply allows to represent more efficiently the full call-graph.

E. Exporters
Along with framework support, the plugin architecture of qprof is also used to implement exporters that will transform the abstract quantum program representation described in Section II D 2 to a more usable format.
Just as the framework plugins, exporter plugins should implement a specific interface schematised in Figure 6. qprof natively implements two textual exporters: one that outputs a gprof-compatible format and another that returns a JSON-formatted string that directly represents a flat call-tree structure used internally by the qprof-compatible exporter.

Flat call-tree representation
Before the profiler report generation, it is convenient to summarise the information contained in the generic callgraph structure presented in Section II D 2. To do so, the gprof and JSON exporters both rely on a flat structure that represents a directed call-tree (i.e. a directed callgraph without loops).
This structure puts an additional restriction to the quantum programs that can be profiled using these exporters: the interdiction to have recursive subroutines (a subroutine that ends up calling itself). It is important to realise that this restriction does not have a huge impact on the area of application of qprof because, as of today, recursive subroutine calls are not widespread in the quantum computing programs and the restriction only applies to the gprof and JSON exporters, the core logic of qprof being capable of handling recursive subroutine calls.
The flat call-tree structure stores, for each subroutine A encountered in the call-graph exploration, a list of all the subroutines B called by A. Along with each called subroutine B, the structure stores the number of times B has been called by A and the execution time that was needed for these calls. Finally, in order to simplify the report generation, each called routine B will also store a list of the routines A it has been called by. Within this list is also stored the number of calls to B that have been performed from each A and the execution time of these calls.

gprof output
The gprof exporter aims at generating a profiler report that is compatible with the profiler report returned by gprof, a well-known classical profiler. Being compatible with a tool that has been around for decades and is still actively used has several advantages.
First and foremost, the fact that a tool that has been stable for decades and is still actively used shows that it provides satisfaction to its users, meaning that the output format includes enough information and is sufficiently easy to read and use in practice.
Secondly, a decades-old, largely used, output format is likely to have a lot of official or user-contributed tools to help analysing and representing it in the best way possible. This is the case for the gprof format that can be translated to a call-graph using the gprof2dot tool and the dot executable from the Graphviz library.
Finally, the gprof output is simple to generate: it is a textual file with a simple and regular format.

III. CODE EXAMPLES AND PRACTICAL APPLICATIONS
This section includes several examples of qprof usage on various quantum circuits ranging from a simple Toffoli gate decomposition in Section III A to more complex algorithm implementations such as Grover's algorithm in Section III B. All these benchmarks are performed on circuits generated using the qiskit framework. An example of benchmarking a quantum implementation of a 1dimensional wave equation solver written using the myqlm framework is finally provided in Section III C.

A. Benchmarking a simple program
One of the most simple quantum program that can be benchmarked is the implementation of a Toffoli gate. Such a benchmark has the benefit of being simple enough to be studied by hand which means that we will be able to verify qprof results by re-computing them.
The decomposition of a Toffoli gate as implemented in the qiskit framework is depicted in Figure 7. A com-plete example using qprof to profile the default Toffoli gate decomposition in qiskit is shown in Listing 1.
The output of qprof, which is here in a gprof format with some slight adaptations as discussed in Section II E, can then be analysed. For the sake of readability and brevity, the full gprof-compatible profiler report will not be included verbatim in this paper and will rather be visualised using the gprof2dot tool that allows representing gprof reports as call-graphs. The call-graph obtained from the report generated in Listing 1 is depicted in Figure 8.
From the call-graph depicted in Figure 8, it is clear that the cost of a Toffoli gate comes from its 6 controlled-X gates, that account for more than 98% of the total execution time. It is also interesting to note that the T gate, known to be very costly when error-correction is needed, is "free" on IBM chips when error-correction is not needed as it is equivalent to a phase change.

B. Grover's algorithm
The Toffoli gate is a good example to start and understand the meaning of qprof's output but the end goal of qprof is to be able to profile large and complex quantum circuits. A good first candidate to show how qprof performs on a more complex circuit is Grover's algorithm.
In this example we use Grover's algorithm on four qubits to find the three quantum states that verify the following formula: The only three 4-qubit quantum states verifying Equation (1) are |0001 , |1001 and |1101 , q 0 being the leftmost qubit in the bra-ket notation.
The code needed to generate the gprof-compatible output for Grover's algorithm with the oracle presented in Equation (1) is given in Listing 2. The resulting call-graph, included in Figure 9, clearly shows that the 100.00% 1× controlled-X gate is still the major contributor to execution time. But this time, contrarily to the Toffoli example shown in Section III A, the controlled-X gate is called by three different subroutines that all contribute significantly to the overall cost: c3z, ccz and mcx. Thanks to qprof it is now easy to understand the subroutines that contribute the most to the total execution time. More importantly, the qprof-compatible report and the call-graph representation gives very insightful information about subroutines calls that are crucial for circuit optimisation. Such information can be used to weight the impact of a given optimisation and then decide whether or not it is worth applying it.
For example, knowing that the ccz subroutine takes 18.61% of the total time, it is easy to deduce that a 20% improvement in the implementation of ccz will translate into a tiny 18.61% 5 = 3.72% improvement to the overall execution time, which might not be worth the effort. On the other hand, optimising the c3z subroutine to reduce its execution time by 20% improves the overall running time by 9.22%, which is nearly 10% and might be an interesting optimisation target. Finally, the call-graph visualisation conveys clearly the information that the cx gate is the most costly subroutine of the Grover's circuit, meaning that even a slight optimisation of the cx execution time will have a high impact on the overall implementation run-time.

C. Quantum wave equation solver
Finally, we include in this paper a more complex example that has been implemented in a previous work with myQLM, a quantum computing framework maintained by Atos. The code used to generate the benchmarked quantum program is available at https://gitlab.com/ cerfacs/qaths/ and is explained in [27].
This example demonstrates that, as can be seen in Figure 10, qprof interface stays nearly the same even though the framework used is now completely different. The only exceptions are some additional parameters (such as linking set in Figure 10) that are directly forwarded to the framework plugin used and additional gate definitions in the gate times data structure because of the way gate decomposition is handled in myQLM.
The call graph obtained by running Figure 10 is reproduced in Figure 11. In order for the call-graph to be readable on a paper format, some subroutines and calls (i.e. nodes and edges respectively) have been discarded from the graphical representation. The call-graph clearly shows that most of the execution time is spent in the oracle implementation. Moreover, (multi-)controlled-X gates are the major contributors to the total execution time. FIG. 9. Call-graph for the Grover's algorithm implementation. Quantum gates included in the gate times variable (here u, u1, u2, u3 and cx) are considered as native gates.
Only native gates have a non-zero self-time as they are the only gates that are really executed on the hardware. Some percentages might not add up to exactly 100% due to rounding errors.

IV. CONCLUSION
In this paper we introduced an open-source and, to the best of our knowledge, novel library that is able to generate profiling reports in well-known formats from a quantum circuit implementation. Our library is able to natively read quantum circuits from multiple frameworks -currently Qiskit and myQLM -and can be easily extended to support more quantum computing libraries. It generates consistent reports independently of the underlying framework used. qprof opens new optimisation opportunities for quantum scientists and programmers by allowing them to view their quantum circuit implementation in a well-known, familiar, synthetic and visual representation. In this paper, we first presented the main concepts used in the internals of qprof and required to make it as extendable as possible: the plugin interface used to wrap different quantum computing frameworks, the call-graph data structure and finally the different native exporters. Then, we used qprof on three different quantum circuit implementations of increasing complexity to demonstrate its features: simplicity of use, adaptability, from qaths.applications.wave_equation.evolve_1D_dirichlet import evolve_1d_dirichlet from qaths.applications.wave_equation.linking_sets.arithmetic_adder import get_linking_set as arith_linking_set from qprof import profile # Circuit generation time = 0.1 discretisation_size = 2 ** 10 epsilon = 1e-3 trotter_order = 1 routine = evolve_1d_dirichlet(time, discretisation_size, epsilon, trotter_order) # Gate execution time definition G = {"u1": 0, "u2": 89, "u3": 178, "cx": 930} gate_times = { "cu1": 2 * G["cx"] + 2 * G["u1"], "cu2": 2 * G["cx"] + 2 * G["u3"], "X": G consistency of the interface and efficiency.
In the future, we plan to extend the set of supported quantum computing frameworks. The number of exporters can also be improved to handle different output formats such as a perf event [12] compatible format, allowing to easily use new visualisations such as flame graphs [18].

SUPPLEMENTARY MATERIAL
The qprof tool is available at https://gitlab.com/ qcomputing/qprof/qprof. FIG. 11. Call-graph of the quantum wave equation solver. Nodes (i.e. quantum routines) that account for less than 0.5% of the total execution time are not plotted. Edges (i.e. subroutine calls) that account for 0.1% or less of the total execution time are also discarded for readability purposes.


Outlier Detection in Smart Grid Communication

Industrial Control System (ICS) networks transmit control and monitoring data in critical environments such as smart grid. Cyber attacks on smart grid communication may cause fatal consequences on energy production, distribution, and eventually the lives of people. Since the attacks can be initiated from both the inside and outside of the network, traditional smart grid security tools like firewalls or Intrusion Detection Systems (IDS), which are typically deployed on the edge of the network, are not able to detect internal threats. For this reason, we also need to analyze behavior of internal ICS communication.

I. INTRODUCTION
ICS communication provides proper functioning of critical infrastructure systems. These systems are naturally exposed to external threats including cyber attacks [1]. Traditionally, industrial systems are well protected against external threats through the use of firewalls and IDS devices that filter communication between ICS systems and Internet traffic, making direct attacks on smart grid communication a rare case. However, attackers can gain access to the system by sending malware to a user via an infected e-mail attachment [2].
The importance of the research is motivated by recent cyber attacks against critical infrastructure systems. One of the attacks against Ukrainian power company happened in December 23, 2015 when BlackEnergy malware caused power disruption to 225,000 customers, lasting up to 6 hours [2]. This attack included multiple stages starting with "spear phishing" e-mails targeting a staff to gain access to the corporate network of the power company. Once inside the power company network, attackers gathered credentials and used VPNs to get access to the internal network. More recent ransomware attack against the Colonial Gas Pipeline in the U.S. happened in May 2021 [3]. Similarly to BlackEnergy attack, it was also initiated from an infected internal station.
Compared to standard information and communication systems, smart grid communication exhibits stable, periodical, and regular communication patterns since the communication occurs between devices with no or little human interference. Typically, a controlling station periodically requests status data from a field device like the Programmable Logical Controller (PLC) or Remote Terminal Unit (RTU) in order to provide a real-time view on industrial processes.
To detect internal threats, we need (i) to regularly monitor smart grid communication and (ii) observe suspicious patterns that occur in the network traffic. One solution is to employ ICS monitoring using extended IPFIX protocol [4] that retrieves monitoring data about active ICS communication. To detect unknown adverse events or unusual behavior, we can observe statistical patterns using ICS flow data.
In this research, we closely examine statistical distribution of inter-arrival times of IEC 104 packets. This is a part of my PhD study that is focused on anomaly detection of ICS communication using statistical methods. The main idea behind the research comes out from stable communication patterns that are typical for ICS communication and can be observed on packet and flow levels. Within the PhD. research, we plan to apply various statistical methods on smart grid communication and evaluate their efficiency for covering common ICS attack vectors defined by MITRE ATT&CK for ICS matrix 1 .

A. Contribution
This research paper presents a proof-of-concept to outlier detection in smart grid communication. We observe packet inter-arrival times of IEC 104 communication and create a statistical profile of the normal traffic. Then we apply the Local Outlier Factor (LOF) algorithm to detect outliers which represent anomalies. Our preliminary results proves viability of this approach for statistical-based anomaly detection in ICS communication.

B. Structure of the Paper
This paper is structured as follows. In Section II we give an overview of published works related to statistical anomaly detection in ICS networks. Next, we briefly introduce IEC 104 protocol, statistical features and the LOF algorithm. Section IV shows our preliminary results. The last section concludes the work and discusses future steps.

II. RELATED WORK
Exploitation of timing attributes such as average packet inter-arrival times and the number of packets or bytes transmitted in a certain interval for anomaly detection has been studied in the past. Barbosa et al. [5], [6] investigated the use of spectral analysis to uncover traffic periodicity. Udd et al. [7] examined a TCP sequence prediction attack for the IEC 104 protocol. Other works mostly focused on flooding and DoS attacks [8]. We explore what attack vectors described by MITRE ATT&CK for ICS can be covered by LOF method.
To enhance detection ability of cyber attack, deep packet inspection (DPI) implemented in IDS tools and timing models have been proposed. Sayegh et al. [9] modelled the inter-arrival times between signatures (i.e., packet sequences) and validated this technique with large amount of injected signatures. Barbosa et al. [10] proposed to model the period of repeated requests in an orderless group. They evaluated the approach on Modbus and MMS datasets without attacks and set relaxed thresholds to avoid high false positive rates. This prevents the detection of subtle changes within a single period which is covered by our approach.
More recently, sequence-aware approaches have been explored. Yang et al. [11], Goldenberg and Wool [12], and Kleinmann and Wool [13] used finite automata to model message sequences of IEC 104, Modbus TCP and S7 respectively. Casselli et al. [14] modelled a sequence of messages using discrete-time Markov chains in order to detect sequence attacks. These approaches observe the order of messages. Since they require deep packet analysis, they have high demands on processing. Modeling of statistical behavior as proposed in this paper is fast and simple for implementation while giving comparable results.

A. Smart Grid Communication
Industrial protocols used in smart grid communication include protocol IEC 61850 (GOOSE, MMS), Modbus, IEC 104, DNP3, DLMS, and others [15]. These protocols transmit control and status data from industrial processes running on RTUs or IEDs. Protocols like IEC 104, DNP3, MMS or DLSM communicate using a client-server model. A master (controlling) station sends commands to a RTU slave (controlled station) in control direction while the slave delivers monitoring data in monitor direction. Protocols like GOOSE or Modbus use a publish-subscribe mechanism in which an application writes data into a local buffer, which is then periodically transmitted to a subscribed agent via L2 multicast.
1) IEC 104: IEC 104 protocol is an application protocol that consists of Application Protocol Control Information (APCI) and Application Service Data Unit (ASDU) sub-layers. It is implemented over TCP but for monitoring purposes we observe so-called virtual flows that represent records with a single ASDU packet transmitted on wire. We focus on interarrival times but there are additional attributes that can be used for statistical model: ASDU size which is also stable for specific IEC 104 commands, frequency of selected commands, e.g., spontaneous events, activations, packet size, etc.
2) Packets inter-arrival time: Packets inter-arrival time is the time taken between the arrival of two subsequent packets. It is computed as a difference between timestamps of two these packets. Inter-arrival time is characteristic for ICS transmissions and exhibits stable and predictable patterns. Changes in inter-arrival times indicate an anomaly on the communication channel. As demonstrated in our previous research [16], additional features like packet size or modeling of exact IEC 104 commands do not improve the accuracy of the statistical approach. We plan to confirm this observation by future experiments. For modelling inter-arrival time of industrial protocols with master-slave communication profile, it is natural to observe statistical distribution of bi-directional transmissions as depicted in Fig. 1.

3) Dataset:
In our experiments, we work with the IEC 104 protocol that is commonly used in smart grids for substation control. We utilise four datasets obtained from the Brno University of Technology testbed, see Tab. I. In order to test the ability of attack detection we created a set of ICS flow records with emulated attack and a AD testing tool 2 .

B. Statistical Outlier Detection
Outlier detection is a statistical technique that discovers data points that are inconsistent with the rest of the data. In general, outliers and inliers are determined using a data distribution model. This paper focuses on unsupervised outlier detection.
Unsupervised methods do not require data labels. As a result, they are more adaptable. The fundamental idea behind unsupervised outlier detection is to score data points solely on the essential characteristics of the dataset. In general, density or distance are used to determine whether a data point is an inlier (normal) or outlier (anomaly). In this proof-of-concept study, we apply Local Outlier Factor (LOF) method [17].
The scientific research on statistical outlier detection provides two approaches to handle outliers in a dataset. First, outliers must be identified for further investigative process.
Second, the data model should be designed to handle outlier data points accurately.

C. Detecting Outliers in Smart Grid
Outlier detection is a statistical procedure that finds suspicious events or items that differ from a dataset's normal form. Outliers are detected as data points that have a significantly lower density than their neighbors. The purpose of outlier detection is to detect rare events or unusual activities that differ from the majority of data points in a dataset [18]. Outlier detection can detect global or local outliers. For a global outlier, outlier detection considers all data points, and the data point is considered an outlier if it is far away from all other data points. The local outlier detection covers a small subset of data points at a time. A local outlier is based on the probability of data point being an outlier as compared to its local neighborhood which is measured by the k-Nearest Neighbors (kNN) algorithm. LOF is a densitybased unsupervised anomaly detection method that computes a given data point's local density deviation with respect to its neighbors. LOF scores are computed for all data points according to parameter k (the number of nearest neighbors) as follows [19]: Definition 1: d(p, o) is the Euclidean distance between two data points p and o. The distance between two data points p and o is calculated using an Euclidean n-dimensional space: Let D be a dataset and k a positive integer. For a data point p, the k-distance (p) is the distance d(p,o) between p and the farthest neighbor data point o by the following conditions: 1) At the least, k data points (records) o ∈ D\{p} maintains that d(p, o ) ≤ d(p, o). 2) At the most, k − 1 data points (records) o ∈ D\{p} maintains that d(p, o ) < d(p, o). Definition 2: k-Nearest Neighbors of p. The meaning of k-Nearest Neighbors of p is any data point q whose distance to the p data point is not greater than the k-distance(p).Those k-nearest neighbors of q form the so called k-distance neighborhood of p, as described below: Definition 3: Reachability distance of p with respect to o. Let k be a positive integer. The reachability distance of data point p with regard to o is as follows: The principle of LOF reachability is depicted in Fig. 2, further details are available at [17].

IV. PRELIMINARY RESULTS
We created a prototype 2 tool implementing LOF algorithm. We applied LOF on IEC 104 inter-arrival times in order to check if LOF modeling is suitable for ICS communication. We used LOF algorithm to learn inter-arrival time distribution of the transmitted IEC 104 packets by computing the k-distance, reachability distance and density of data points. The algorithm uses the learned model to raise an alarm for packets that significantly differ from the model. Figures 3 and 4 present a graphical representation of inter-arrival times of two datasets. The red points denote inter-arrival time distribution detected as anomaly (outliers) while the blue points represent normal communication. We discovered that applying LOF to a large data stream is extremely computationally inefficient and may lead to incorrect prediction results. We then applied the LOF on time windows where each data block contains 5000 inter-arrival time records.
The main advantage of the LOF algorithm is that it works well with stream data and detects outliers with respect to Fig. 4. Anomaly detection in mega104-14-12-18 PCAP file density of their neighbouring data records. Also, this algorithm is able to detect outliers regardless of the data distribution of normal behavior, since it does not make any assumptions about the distribution of data records. Our preliminary results show that the proposed method works well with ICS data and is easy to operate.
Validation tests were applied to determine the viability of the proposed technique. We tested the stability of the approach using our datasets. For our validation tests, we divided normal communication into two parts. We used the first two thirds of the communication to train the model. Then, we tested these ranges on the last third of data. After testing, we classified all testing packets as normal communication. The results demonstrate high accuracy of the proposed approach. The method did not produce any false positive during the test.

V. CONCLUSION
In this paper, we presented a proof-of-concept method for anomaly detection of smart grid control protocols. We observed inter-arrival times of IEC 104 communication and applied LOF algorithm for outlier detection. The preliminary results show that LOF creates a stable statistical model for ICS traffic and is able to detect outliers caused by cyber attacks. The proposed security mechanism was tested and validated using our IEC 104 datasets. Experiments show the suitability, usability and high accuracy of the proposed method on smart grid communication. In the future work we plan to apply this approach to other ICS protocols like GOOSE, MMS and DLMS. We also plan to observe other statistical features like packet or flow size. Although our key goal was to test the applicability of the proposed solution, we also tested the performance of the developed tool. We gathered basic time statistics. The processing time of extracting key attributes from the PCAP file into a CSV file took an average of one minute and 10 seconds. This time had small deviations because it takes time to initialize Tshark. The time to learn a CSV file depends on the number of packets. At average, to learn 5,000 packets, 20 seconds were required.
The results are part of my PhD research that focused on statistical-based anomaly detection in ICS protocols. In the future work, we focus is extending detection model to common cyber attacks defined by MITRE ATT&CK ICS matrix.

ACKNOWLEDGMENT
The work is supported by the Brno University of Technology project "Application of AI methods to cyber security and control systems", no. FIT-S-20-6293.


extfieldsofstudy:
Computer Science
Physics

text:
Performance Evaluation of Digital Transmission with Matlab on Eye Pattern Parameter

The effect of digital transmission impairments is seen from a diagram called an eye pattern. Measurements of the device can be represented in a simulation. Our motivation in conducting this simulation aims to dig deeper into the variation of parameter changes that exist in improving the eye pattern opening angle. The raised cosine pulse with roll-off (a) of 0.3 was used as the basis for this study, using the Matlab simulation. Important parameter values have been explored to show the best data. Finally, we suggest that we continue to explore optimization techniques for the best results for this eye pattern visualization.

Introduction
Distortion during transmission of telecommunication is a fact of life. One source of distortion to a pulse train is random noise, and we will take up this in the next section. A host of non-noise distortions, or "impairments" as they are called, can affect signals. Channel filters are a necessity because they make a signal bandwidth-efficient and reduce interference to neighboring transmissions, but filters also distort signals. So also do nonlinearities in signal amplifiers. Another kind of distortion occurs when the timing of the sampler is faulty so that the signal there is sampled too early or too late. Some of these impairments are supported with the aid of a convenient tool called the eye pattern. [1] Eye patterns are important in the operation of digital transmission devices because they can help monitor device performance or as a parameter checkpoint when troubleshooting digital transmission devices. In addition to eye patterns, the constellation pattern modulation process includes this checkpoint parameter. [2] The aim of this research is to simulate and visualize eye patterns as an important parameter of digital transmission equipment. The simulation is carried out by comparing the results of various changes in device system parameters so that an analysis of the opportunities for optimizing device performance can be carried out on existing parameters. With the help of Matlab on this eye pattern and related parameters, it is possible to deepen monitoring of the performance of digital transmission equipment in the future.

Raised Cosine Pulse
Much of what happens to pulse trains can be characterized as inter-symbol interference [ISI]. Loosely defined, ISI is the effect of a pulse on the detection in other symbol intervals. Figure 1 shows the standard 30% excess-band width raised cosine (RC) pulse together with a delayed version of the same pulse and a low-pass filtered version. The original pulse is Nyquist. Consequently, in a sampling receiver with accurate sampling times nT, the pulse is zero at all sample times other than its own and does not affect the other sampler outputs. In the linear receiver, with root RC pulses, the response at the receiver filter output to each root RC pulse is the RC pulse, and here, again, with accurate sample timing, each transmitted pulse affects only one sampler output. The second pulse in Figure 1 is delayed by 0.2T so that the sampling occurs 0.2T too early. Now the pulse contributes at other sampling times, in the amounts shown in the figure. This is ISI stemming from sampling time error. The third pulse is the RC pulse, filtered at a cut-off frequency of 0.5 / T Hz by a six-pole Butterworth filter. The filter has a delay of around 1.4 T seconds; but even with this taken off, at the other sampling times, substantial contributions exist. This ISI stems from the time dispersion of the filter.

Eye Pattern
A classical technique for visualizing the effect of Inter Symbol Interference (ISI) is the eye pattern. The eye pattern gives a qualitative measure of system performance. A well-defined and open eye usually indicates good performance, while a poorly defined eye usually indicates poor performance. In addition, the size of the eye relates to the accuracy required of the symbol synchronizer. While the eye diagram does not provide a quantitative measure of system performance, it is difficult to conceive of a high-performance system having a poorly defined eye diagram. Three segments of a waveform, with each segment corresponding to a symbol period, are shown in figure 2. The waveform corresponding to three data symbols is illustrated in figure 2(a). Assume that this waveform is displayed on an oscilloscope and that the oscilloscope is triggered at the points denoted by the dotted vertical lines. The result will be the three-segment eye diagram illustrated in figure 2 Figure 2. Generation of an eye pattern. [3] ISI is constructed by overlapping multiple segments of the received waveform over a fixed window, which tells us how different combinations of symbols could potentially create ISI. For an ideal channel and square root Nyquist pulses at either end, the eye is open. However, for the dispersive channel, the eye is closed. An open eye implies that, by an appropriate choice of sampling times, we can make reliable single-sample symbol decisions, while a closed eye means that more sophisticated equalization techniques are needed for symbol recovery. [4] The eye pattern is obtained by superimposing a large number of trajectories, with duration 2T, of the matched filter's output signal. It can be displayed on an oscilloscope by synchronizing with the symbol rate 1/T. The afterglow of the screen makes it possible for the superposition to persist. In the presence of noise, as we are going to see, the wider the eye is vertical, the lower the error probability. Therefore, we have to choose the decision time where the eye is vertically "widest". The low level of interference is clearly shown by the eye pattern without noise represented in figure 3. The trajectories almost converge to the same point at multiplies of T. Therefore, if the sampling is done at these times, the values located around 1 are likely to correspond to the transmission of a bit 1. Eye pattern data is a representation of a high-speed digital signal that enables quick visualization and determination of key parameters of the electrical quality of a signal. [  A method use least-median-of-squares (LMS) location estimator for analyzing eye patterns always provides a unique solution stated in [8]. The LMS protocol is indifferent to outliers and data distributions, as opposed to widely used histogram approaches. The reason for designing this algorithm is to create an autonomous, benchmark approach that is both suitable for detailed study of ambiguity and can act as a comparative tool as there are currently no standardized industry algorithms. Using this method, the basic parameters of an eye diagram are determined by the authors, namely the levels of one and zero, and the time and amplitude crossings. With these parameters determined, various efficiency measures, such as extinction ratio and root-mean-square jitter, can be obtained and the eye-mask synchronization performed.

Method
In this study, the changes in sequence length of the source and binary velocity were observed to the values of other parameters physically by observing the existing eye patterns. The tool used is Matlab. Although Matlab has its own eye diagram routine, a code snippet is provided here to convey the concept of this research clearly. The output from the resulting receiver filter is the input to these fragments, but in general, it can be plotted an eye diagram based on the baseband waveform at each stage in the system. Here the two roll-off prices α are taken as 0.3 and compared with 0.6 and 1.0 at several SNR values. The modulation scheme used is QPSK. Previously, in the baseband section, observations were made on level four PAM signal data. The complete parameters that the authors take in this study are listed in table 1.

Result and Discussions
In this research, a PAM signal with two bits per symbol was used. Therefore obtained the alphabet symbol set is (-3 -1 1 3) as figure 7. The overall impulse ensures that the ISI is low enough for the symbol-by-symbol decision to yield good results. The increase in SNR taken in the simulation from the 2 dB value which continues to be enlarged, will cause the eye pattern opening to improve, even towards a null noise eye pattern. The constellation of the QAM signal that is formed is shown in Figure 8. At a price of α = 0.3 with a band of the channel of 650 Hz, the eye pattern opening angle still looks small at the SNR price of 15 dB. This is shown in Figure 4, where the qualitative value is similar for the eye pattern price for the price n = 2; B = 800 Hz; α = 0.6; and SNR = 10 dB. In the view of the opening towards the ideal, it can be seen in Figure 5. This can happen at other prices such as at B = 900; α = 0.8; SNR = 30 dB. The other side of the increase in eye pattern openings occurs with an increase in the price of B which also results in an increase in α. For the simulation, the increase in α corresponds to the increase in the band of the channel. Meanwhile, the increase in the angle of the eye pattern is getting bigger as the SNR is taken. Conversely, when the binary rate is reduced to

Conclusion
The effect of the eye pattern opening depends on the amount of SNR in the system, the bigger the SNR, the cleaner the eye pattern. In line with that, an enlarged α value at a fixed SNR will also produce a better eye pattern opening.