Ioannis Kyriakidis
email: kyriakidis@teicrete.gr

Kostas Karatzas

George Papadourakis

Andrew Ware
email: jaware@glam.ac.uk

Jaakko Kukkonen
email: jaakko.kukkonen@fmi.fi

Investigation and Forecasting of the Common Air Quality Index in Thessaloniki, Greece

Keywords: air pollution, common air quality index, neural networks, decision trees, linear regression

Air pollution can affect health and well-being of people and ecosystems. Due to the health risk posed for sensitive population groups, it is important to provide with hourly and daily forecasts of air pollution. One way to assess air pollution is to make use of the Common Air Quality Index (CAQI) of the European Environment Agency (EEA). In this paper we employ a number of Computational Intelligence algorithms to study the forecasting of the hourly and daily CAQI. These algorithms include artificial neural networks, decision trees and regression models combined with different datasets. The results provide with a satisfactory CAQI forecasting performance that may be the basis of an operational forecasting system.

Introduction

High levels of air pollution are associated with a number of health problems and affect quality of life [START_REF]Calculate the CAQI Levels from forecasted Daily/Hourly CAQI values by using regression models[END_REF]. As an example it may be noted that ambient air pollution, in terms of fine particulate matter (PM2.5), causes about 3% of mortality from cardiopulmonary disease, about 5% of mortality from cancer of the trachea, bronchus, and lung, and about 1% of mortality from acute respiratory infections in children (age<5yr), worldwide [2]. On this basis, the importance of forecasting of poor air quality is evident.

Two forecasting approaches may be employed. The first one makes use of deterministic, numerical models that are capable of "representing" the physical and chemical processes affecting air pollution. They have the advantage of fully covering a geographical area of interest, but they require detailed input in terms of emission data, meteorology and physico-chemical parameters, and they are computationally demanding. The other forecasting approach makes use of data-driven modeling, employing either statistical [START_REF] Vlachogianni | Evaluation of a multiple regression model for the forecasting of the concentrations of NOx and PM10 in Athens and Helsinki[END_REF] or Computational Intelligence (CI) methods [4], [5].

We make use of the CAQI of the EEA as the parameter of interest that needs to be forecasted in order to provide health related warnings and information to people. We selected the city of Thessaloniki, Greece as our target city, and we made an effort to construct operational forecasting models, with the aid of CI methods that include Artificial Neural Networks (ANNs) and Decision Trees (DT), as well as statistical methods (Linear Regression Models). In the rest of the paper we are describing the materials and methods that we used (chapter 2), the results of our calculations and the relevant discussion (chapter 3), while we draw our conclusions in chapter 4.

Materials and Methods

Area of interest and datasets used

Thessaloniki is the second largest city of Greece, and is characterized by a pronounced problem of air pollution, especially in terms of Particulate Matter (PM) and more specifically of respirable particles (PM 10 , i.e. particles of aerodynamic diameter smaller than 10 μm). Air quality monitoring is conducted in various locations of the Greater Thessaloniki Area (www.airthess.gr), and data from these measurements were used in the frame of the current study.

We analyzed data from four monitoring stations located in (a) Agia Sofia (city center), (b) Panorama (a suburb in the northern mountainous area of the city), (c) Sindos (industrial-suburban area in the west of the city influenced by traffic) and (d) Kordelio (a densely populated area close to the industrial area of the city, influenced by traffic). These data included hourly measurements of air pollutants such as Carbon Monoxide (CO), Nitrogen Dioxide (NO 2 ), Nitrogen Monoxide (NO), Ozone (O 3 ), Respirable particles (PM 10 ), and Sulphur dioxide (SO 2 ) and they were complemented by meteorological observations. Data came from the years 2001-2003 (time period with a low percentage of problematic or missing data).

The Common Air Quality Index

There has been an initiative to combine the impacts of various pollutants and come up with an aggregated measure that may be used for the description of the atmospheric quality. Various air quality indices have been suggested by institutes and authorities in various countries [6], [7]. The European Environment Agency has adopted an index called the Common Air Quality Index (CAQI, www.eyeonearth.org) to identify the atmospheric quality towards public safety and health.

The CAQI is designed to present and compare air quality in near-real time on an hourly or daily basis. The CAQI has 5 categorical levels, corresponding to a range of values starting from 0 (very low) to >100 (very high). Two types of a CAQI are specified, an urban background index and a traffic related index. The urban background index takes into account three main pollutants (NO 2 , PM 10 and O 3 ) and two auxiliary pollutants (SO 2 and CO) while the traffic index takes to into account two main pollutants (NO 2 and PM 10 ) and also one auxiliary (CO). Those are described in detail in [8], and they are calculated as explained in Table 1. This CAQI was used to calculate the air quality levels (hourly and daily) for each of the four locations of interest. The traffic and urban background CAQI were calculated per station in accordance to the nature of the monitoring sites: the former was applied for Agia Sofia, Kordelio and Sindos, while the latter for the Panorama station (Table 2) Table 1. Calculating the two types of CAQI (background and traffic index) by EEA [8]. All concentrations are in μg/m 3 . [CO]: 8hours moving average values

            500 ] [ 5 / ] [ 500 ] [ 100 5 . 37 ) 8 / ] ([ 100 ] [ 0 2 / ] [ ]) ([ _ 2 2 2 2 2 2 2 2 SO SO SO SO SO SO SO Index SO                  400 ] [ 4 / ] [ 400 ] [ 200 50 ) 8 / ] ([ 200 ] [ 100 25 ) 4 / ] ([ 100 ] [ 0 2 / ] [ ]) ([ _ 2 2 2 2 2 2 2 2 2 2 NO NO NO NO NO NO NO NO NO Index NO 5)/12 * ] [O ( ]) _Index([O O 3 3 3  ] [PM ]) _Index([PM PM 10 10 10                   20000 ] [ 200 / ] [ 20000 ] [ 10000 50 ) 400 / ] ([ 10000 ] [ 5000 25 ) 100 / ] ([ 5000 ] [ 0 200 / ] [ ]) ([ _ CO CO CO CO CO CO CO CO CO Index CO Main Background Index = Max(NO2_Index, PM10_Index, O3_Index) Auxiliary Background Index = Max(NO2_Index, PM10_Index, O3_Index, SO2_Index, CO_Index) Main Traffic Index = Max(NO2_Index, PM10_Index) Auxiliary Traffic Index = Max(NO2_Index, PM10_Index, CO_Index)

Pre-processing of input data and methods for the forecasting of the CAQI

The categorization of the CAQI into various levels reveals that in order to forecast the CAQI for the next hours and days, we need to make use of methods that can effectively classify existing (historical) data, and forecast the classification value of interest. On this basis, and taking into account relevant literature, we employed CI methods, and more specifically Artificial Neural Networks and Decision Trees. ANNs were advanced in the late '80s, popularizing techniques like Multi-layer Perceptrons (MLP) [9] and Self-Organizing Maps (SOM) [10] while they can be trained to successfully approximate virtually any continuous function [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]. Their advantages also include greater fault tolerance, robustness, and adaptability, especially compared to expert systems, due to the large number of interconnected processing elements that can be trained to learn from new patterns [START_REF] Lippman | An introduction to computing with neural nets[END_REF]. These features provide ANNs with the potential to model complex non-linear phenomenon like air pollution [START_REF] Chelani | Prediction of ambient PM 10 and toxic, metals using artificial neural networks[END_REF].

A Decision Tree on the other hand, is a hierarchical data structure implementing the divide-and-conquer strategy [START_REF] Ethem | Introduction to machine learning[END_REF]. It is an efficient nonparametric method, which can be used for both classification and regression. DTs are essentially a map of the reasoning process in which a tree-like graph is constructed to explore options and investigate the possible outcomes of choosing the options. The reasoning process starts from a root node, transverses along the branches tagged with decision nodes, and terminates in a leaf node [START_REF] Lakhmi | Advances in Intelligent Methodologies and Techniques[END_REF].

The above two CI methods were used together with Linear Regression (LR) models, which are based on a linear combination of the input parameters. Such models may take into account the "memory" and inertia interwoven to the air quality system of any area and have been proven successful especially if they are "fed" with lagged values, i.e. past values of the parameter of interest.

From the results of Table 2 it was made evident that Agia Sofia station demonstrates the highest values in terms of "High" and "Very High" CAQI classes. For this reason, it was selected as the station for the development of CAQI forecasting models. Air quality and meteorological data used for LR, ANNs and DT model construction, originate from this location. In order to forecast hourly and daily CAQI, three datasets were used, depending on the forecasting target. daily, according to the type of index being calculated. In addition, this dataset also included air quality and meteorological values. Dataset 3 consists of 1042 daily records and 23671 hourly records. Concentration data were preprocessed before used as input to the forecasting algorithms (ANNs and DTs). The first step of the preprocessing procedure was to remove all records that did not include values for the target characteristic(s). Specifically, for Datasets 1 and 3, all records that do did include CAQI values, were removed. In the case of Dataset 2 all records that did not include all target values (CO, SO 2 and O 3 ) were also removed. The next step was to use an interpolation method to calculate all missing values [START_REF] Kyriakidis | Using Preprocessing Techniques in Air Quality forecasting with Artificial Neural Networks[END_REF]. The wind direction values were encoded by using the formula WD=1+sin(θ+π/4), in order to replace the cyclic nature of this variable with a linear one. In addition, as ANNs require numerical values as input, nominal values were converted to numerical values for the CAQI, as described in Table 3. The specific numerical values were selected, because of the hyperbolic tangent sigmoid transfer function that was used to normalize the input data in the case of ANNs and DTs. Table 4 presents the way that the continuous numerical values (forecasts) were mapped back to nominal CAQI values. For this mapping the numerical values were divided in five equal in range groups. Very High

The performance of the models was estimated with the aid of statistical indicators like the Index of Agreement and the Cohen's Kappa, which are defined as follows:

            n i i i n i i i O O O P O P IA 1 2 2 1 1 (1)
where P i and O i are the forecasted and observed concentration values respectively. The index of agreement ranges between 0 and 1, with 0 stating no agreement, and best agreement given by the limiting case of 1.

The Cohen's Kappa or critical success index (CSI) verifies how well the events of interest (CAQI values) were forecasted and it is unaffected by the number of correct forecasts. Its calculation is based on the formulae:

c b a a CSI   
where a,b, c and d are calculated as follows: a= true positives, b= false negatives, c= false positives, d=true negatives. Positives or negatives refer to the occurrence or not of the event of interest (i.e correct classification of the CAQI), respectively. The CSI used in this study does not use the "d" events (true negatives). As a consequence, the CSI is tending to give poorer scores for rare events, but it takes into account both false negative events and false positive events. Therefore, the CSI can be characterized as a more balanced score [START_REF] Ji | Statistical comparison of interplanetary conditions causing intense geomagnetic storms (Dst<=-100 nT)[END_REF].

Results and discussion

CAQI Forecasting with the aid of Linear Regression Models

Forecasting of the numerical hourly and daily CAQI values

LR models were used as the reference method for the development of the CAQI forecasting models. For this reason lagged values of the CAQI (a total of up to 10 values, Dataset 1) were tested. Results indicate that the hourly CAQI numerical values for Agia Sofia can be forecasted by a LR model by employing only one lagged value (IA = 0.92). Contrary to the hourly, the daily CAQI numerical values were not as successfully forecasted (IA = 0.76).

Forecasting of the nominal hourly and daily CAQI levels

The calculations were repeated for the CAQI levels (nominal values), in addition to their numerical values. Overall results are presented in Tables 5 and6. On this basis it is made evident that both hourly and daily CAQI levels are not forecasted in a satisfactory manner with the aid of LR models. This can be verified via their low Cohen's Kappa Index (Table 7). As a next step, the performance of the LR models for the next day's CAQI levels forecast was further investigated (nominal values). For this purpose, the equal as well as the weighted participation of hourly CAQI values in the forecasting model was employed, as described below.

Forecasting of the numerical 24-hourly values to calculate the daily CAQI levels

The fact that the forecast of hourly CAQI values reached a 0.92 IA led us to investigate the forecasting of the daily CAQI levels via forecasting first the 24 hourly values of the next day. To do so, the hourly CAQI values of the previous day were used in the following manner: for the forecasting of the first hourly CAQI value for the next day, the 24 hourly values of the previous day (observed ones) were used. For the forecasting of the second hourly CAQI value of the next day, the previous (forecasted) value and the 23 hourly values of the previous day (observed ones) were used, and so forth. Via this procedure, the hourly CAQI values were forecasted and then we calculated the daily average CAQI value. Table 8 presents the comparison of daily CAQI forecasting performance via LR models at the station of Agia Sofia, by using lagged observed numerical values and by (equally) using 24 hourly forecasted values. The hourly CAQI numerical values for the Agia Sogia station achieved an IA = 0.70, on the basis of LR models. The performance was decreased in comparison to the performance of the model that uses ten observed lagged values, which reached IA = 0.92.On the other hand, the daily CAQI Levels performance was increased when using 24 forecasted lagged values with an IA=0.5078. This difference in the performance may be attributed to the use of a higher percentage of observations for the forecasting of the hourly CAQI values for the first hours of the next day, thus allowing information of the previous day to "penetrate" the target day of the forecast.

Daily CAQI Levels based on 24-h Forecasted Values with Weighting Factors

In order to further investigate the influence of input values to the performance of the forecasting models, additional daily CAQIs were calculated by using weighted factors. These weighted factors dictated that the first forecasted values of each day (values at midnight, 1 a.m., 2 a.m. etc.) will contribute more than the last forecasted values. Table 9 presents one of the weighting methods tested, named "Factor4", which provides the best forecasting performance, in comparison to others. When using the hourly forecasted CAQI values to calculate the daily CAQI, the results are slightly better (IA= 0.5105, i.e. an improvement of 1,51%) than those when forecasting directly the daily CAQI.

CAQI Forecasting with the aid of ANNs and Decision Trees

In order to develop ANN and DT models, data were firstly normalized by applying the hyperbolic tangent sigmoid transfer function, which was also applied on the hidden layers. On the output layer the linear transfer function was used. This is a common structure for function approximation (or regression) problems. The first 50% of the data is used for training, the next 25% for validation and finally the last 25% for testing. Validation data are used in order to early terminate training (i.e. if the network performance on the validation data fails to improve), thus avoiding overfitting. Test data are used for the estimation of the model performance and the calculation of the relevant indices. Concerning the ANN architecture, one hidden layer was used in all models. The number of neurons for each model was equal to the number of the input parameters. For the training phase we implemented the Levenberg-Marquardt Back propagation algorithm [START_REF] Saini | Artificial neural network based peak load forecasting using Levenberg-Marquardt and quasi-Newton methods[END_REF]. Principal Component Analysis (PCA) [START_REF] Voukantsis | Intercomparison of air quality data using principal component analysis, and forecasting of PM 10 and PM 2.5 concentrations using artificial neural networks[END_REF] was performed to the data in order to change their dimensionality, as this helps the ANN to find patterns into the data. It should also be noted that none of the input parameters were excluded. Several experiments were performed in order to forecast the CAQI levels by using ANNs and Decision Trees:

1 .

 1 Dataset 1: Dataset with only lagged (i.e. previous) CAQI values (from 1 to 10), either hourly or daily, according to the type of index being calculated. Dataset 1 consists of 1042 daily records and 23671 hourly records. 2. Dataset 2: this included air quality (CO, SO 2 , O 3 ) and meteorological values (Temperature, Dew Point, Humidity, Sea Level Pressure, Wind Direction). Dataset 2 consists of 1043 daily records and 23672 hourly records. 3. Dataset 3: this included lagged (previous) CAQI values, (from 1 to 5), hourly or

Table 2 .

 2 The distribution of the CAQI per level and station. The calculation grid for the CAQI is: Very Low (0-25), Low (26-50), Medium (51-75), High (76-100), Very high (>100).

	Index Level	Agia Sofia (2001-2003) Hourly Daily	Panorama (2001-2002) Hourly Daily	Sindos (2001) Hourly	Kordelio (2001-2002) Daily Hourly Daily
	Very high	3510	127	243	1	591	0	2538	95
	High	3706	168	613	12	618	0	1958	112
	Medium	7414	438	5569	306	1917	3	3866	225
	Low	8264	304	9427	377	3558	361	5056	183
	Very low	1666	18	885	21	1309	1	1514	17

Table 3 .

 3 Conversion table from CAQI Levels (nominal values) to numerical values.

Table 4 .

 4 Conversion

	table from continuous
	numerical values to CAQI Levels (nominal
	values).

Table 5 .

 5 Observed vs forecasted hourly CAQI Levels for Agia Sofia via LR. Rows sum 100%.

		Hourly CAQI		forecasted CAQI Levels
			Levels	Very low	Low	Medium	High	Very High
	Observed	CAQI Levels	Very low Low Medium High Very High	65.31 7.14 0.23 0.17 0.03	31.88 74.74 23.22 4.11 1.30	1.10 16.49 61.78 33.55 6.74	0.72 1.18 12.43 46.95 23.00	0.99 0.44 2.34 15.21 68.93

Table 6 .

 6 Observed vs forecasted daily CAQI Levels for Agia Sofia via LR. Rows sum 100%.

	Daily CAQI Levels	Very low	Forecasted CAQI Levels Low Medium High	Very High
	Observed	CAQI Levels	Very low Low Medium High Very High	23.53 5.75 0.64 0.00 0.00	35.29 62.94 33.26 12.87 3.94	41.18 23.64 52.36 47.95 24.41	0.00 6.07 9.66 25.73 29.92	0.00 1.60 4.08 13.45 41.73

Table 7 .

 7 Results of hourly and daily CAQI level forecasts via LR models (Ag. Sofia).

Table 8 .

 8 Comparison of daily CAQI forecasting performance via LR models (Ag. Sofia).

	Table 9. "Factor4"
	weighting of hourly
	values at Agia Sofia.

Conclusions and further research directions

In the present paper the CAQI was used to calculate the pollution levels (hourly and daily) for four locations in Thessaloniki, Greece. On the basis of the results, the Agia Sofia station was chosen as the one for which forecasting models were developed and tested. LR models, ANNS and DTs were used in order to perform the forecasting of hourly and daily CAQI levels (i.e. nominal) and values (i.e. numerical). Different input datasets were used for that purpose, with a different number of lagged values. We observed that when the forecasted daily CAQI value was calculated on the basis of the forecasted hourly CAQI values, the performance was better in comparison to the one achieved via the direct forecast of the daily CAQI values. Moreover, if weighting factors were used in order to calculate the daily CAQI values, the forecasting performance improved. The best performance in terms of the Cohen's kappa was 0.53 for the hourly CAQI and 0.44 for the daily index levels, thus acceptable for operational purposes [19].

to calculate the hourly and daily CAQI (nominal) levels. Table 10 presents the forecasting results of different datasets for hourly and daily CAQI levels, indicated as Model 3 (when Dataset 1 was used), Model 4 (when Dataset 2 was used) and Model 5 (when Dataset 3 was used). 2. Hourly and daily NO 2 and PM 10 values were forecasted with ANNs, in order to calculate the hourly and daily CAQI Levels by using the formulas of Table 1. As input to the ANN, Dataset 2 was used (Model 6 in Table 10). 3. With the aid of ANNs the hourly and daily CAQI Levels were forecasted. In order to convert these numerical values to CAQI Levels (nominal values), the logical expressions of Table 4 were used. Table 10 present the forecasting results of different datasets for hourly and daily CAQI levels, indicated as Model 7 (when Dataset 1 was used), Model 8 (when Dataset 2 was used) and Model 9 (when Dataset 3 was used). 4. With the aid of DTs the hourly and daily CAQI Levels were forecasted. In Table 10 these are indicated as Model 10 (when Dataset 1 was used), Model 11 when Dataset 2 was used) and Model 12 (when Dataset 3 was used).

Comparison of results

On the basis of the overall results presented in Table 10 we can evaluate the forecasting performance for each daily and hourly model was used. For this purpose, we used the "Percentage Agreement" (PA), in order to calculate the percentage of the times that the forecasted CAQI Levels were identical with the observed CAQI Levels, in order to evaluate the models. The best forecasting performance of daily CAQI Levels (PA=61%), was achieved when Model 9 (Forecast Directly the CAQI Levels with Dataset 3) was used with the aid of ANNs, with an increase of forecasting performance reaching 10% compared to the best forecasting performance when LR models were used, and with an increase of 4% compared to the best forecasting performance when DTs were used. The best forecasting performance of hourly CAQI Levels (PA=65%), was achieved when LR models were used, with an increase of 5% compared to the best forecasting performance when ANNs were used, and with an increase of 2% compared to the best forecasting performance when DTs were used.


Ilia Nouretdinov

Dmitry Devetyarov

Brian Burford

Stephane Camuzeaux

Aleksandra Gentry-Maharaj

Ali Tiss

Celia Smith

Zhiyuan Luo

Alexey Chervonenkis

Rachel Hallett

Volodya Vovk

Mike Waterfield

Rainer Cramer

John F Timms

Ian Jacobs

Usha Menon

Alex Gammerman

Multiprobabilistic Venn Predictors with Logistic Regression

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Prediction of heart disease (HD) and ovarian cancer (OC) is a critical task. For some of these diseases (e.g., OC) it is especially crucial in their early stages, when the disease has no clinical symptoms. Mass spectrometry techniques are widely deployed in these problems.

When predicting diagnosis based on proteomics data, very often, the classical machine learning approach is to predict the diagnosis without any measure of how accurate this prediction is. In this work we describe the methodology of proteomics mass spectrometry data analysis based on hedging predictions that is how strongly we believe in this prediction [START_REF] Gammerman | Hedging predictions in machine learning[END_REF].

The framework of Venn machines was introduced in [START_REF] Vovk | Algorithmic learning in a random world[END_REF] and represents a new generation of prediction algorithms. These methods have a range of advantages over the known techniques. Firstly, the prediction which is made is always tailored to the object; as a result, we output a probability interval to each patient's diagnosis. Secondly, the only statistical assumption which is used is the exchangeability assumption which can be satisfied when the data sets are in random order.

Strict definitions of Venn machines are given in Section 2. The main idea is as follows. We first divide examples into categories, the category assigned to an example may depend not only on the example itself, but also on its relation to the rest of examples. For each hypothesis about the new label, we classify the new object into one of the categories, and then use frequencies of labels in the chosen category as predictable distribution of the new object's label. Due to different hypotheses, the machine outputs several (two in the binary case) probability distributions (multiprobability distribution) for the new object rather than the single one.

Practically any known machine learning algorithm can be used as an underlying algorithm in this framework (such as Neural Network in [START_REF] Papadopoulos | Reliable Probabilistic Prediction for Medical Decision Support[END_REF] and SVM in [START_REF] Zhou | A Comparison of Venn Machine with Platt's Method in Probabilistic Outputs[END_REF]). In this work we used logistic regression as an underlying algorithm. It is popular as a method that initially outputs probabilities and provides information about relative weight of features. We will compare Venn machines predictions with probability predictions output by logistic regression.

The methodology is designed for the analysis of proteomic mass spectrometry data collected in the UKCTOCS trial (for more information see www.ukctocs.org.uk).

Venn Machines

Consider a training set consisting of object, x i , label, y i , pairs: (x 1 , y 1 ) , . . . , (x n , y n ). To predict a label y n+1 for a new object x n+1 = x new , we check different hypotheses y n+1 = y each time including the pair (x n+1 , y n+1 ) = (x new , y) into the set.

The idea of Venn machines is based on a taxonomy function A n , n ∈ N, which classifies the relation between an example and the set of the other examples:

τ i = A n+1 ((x i , y i ), {(x 1 , y 1 ), . . . , (x i-1 , y i-1
), (x i+1 , y i+1 ), . . . , (x n+1 , y n+1 )}).

Values τ i are called categories and are taken from a finite set T = {τ 1 , τ 2 , . . . , τ k }. Equivalently, a taxonomy function assigns to each example (x i , y i ) its category τ i , or in other words grouping all examples to a finite set of categories. This grouping should not depend on the order of examples within a sequence.

The conventional way of using Venn ideas was as follows. Categories are formed using only the training set. For each non-empty category τ the following values are calculated: N τ -the total number of examples from the training set, assigned to category τ , and N τ (y ) -the number of examples within category τ having label y . Then empirical probabilities of an object within category τ to have a label y are found as

P τ (y ) = N τ (y )/N τ . (1)
Now, given a new object x n+1 with unknown label y n+1 , one should assign it somehow to the most likely category of those already found using only the training set; let it be τ * . Then the empirical probabilities P τ * (y ) are considered as probabilities of the object x n+1 to have a label y . The idea of conformal predictors [START_REF] Vovk | Algorithmic learning in a random world[END_REF] allows us to construct several probability distributions (multi probability distribution) of a label y for a new object. First we consider a hypothesis that the label y n+1 of a new object x n+1 is equal to y, (y n+1 = y). Then we add the pair (x n+1 , y) to the training set and apply to this extended sequence the taxonomy function A. This groups all the elements of the sequence to categories. Let τ * (x n+1 , y) be the category containing the pair (x n+1 , y). Now for this category we calculate, as previously, the values N τ * , N τ * (y ) and empirical probability distribution

P τ * (xn+1,y) (y ) = N τ * (y )/N τ * . (2)
This distribution depends implicitly on the object x n+1 and its hypothetical label y. Trying all possible hypothesises of the label y n+1 being equal to y, we obtain a set of distributions P y (y ) = P τ * (xn+1,y) (y ) for all possible labels y. These distributions in general will be different, as when changing the value of y we change (in general) grouping into categories, the category τ * (x n+1 , y), containing the pair (x n+1 , y), the numbers N τ * and N τ * (y ). So we obtain, as the output of Venn predictors, as many distributions as the number of possible labels.

In the two-class problem (Y = {0, 1}), Venn predictors have two probability distributions, defined by p y (1) = P {y n+1 = 1}. Thus, the output can be interpreted as the interval

[p - new , p + new ] = [min{p 0 (1), p 1 (1)}, max{p 0 (1), p 1 (1)}] , (3)
which is an estimation of probability that y n+1 = 1. We will refer to p - new and p + new as lower Venn prediction and upper Venn prediction, respectively. They can be interpreted as lower and upper bounds for the probability. Thus if one sets a risk threshold θ and takes all the predictions with lower Venn prediction not smaller than θ then the expected percentage of cases between these examples should be between θ and 1 as well.

A Venn predictor is entirely defined by its taxonomy. In the next section we describe a taxonomy based on the logistic regression.

Logistic Regression

Logistic regression outputs the probability distribution of a new label. It produces these distribution as follows.

Suppose each object out of the training set x 1 , . . . , x n is an m-dimensional vector, each with corresponding labels y 1 , . . . , y n , ∈ Y = {0, 1}.

The statistical model of logistic regression is based on the assumption that P {y i = 1} = 1/ 1 + e -<xi,b> . The optimization goal for logistic regression is:

n i=1 log 1 + e (-1) y i xi,b + a b, b → min b . ( 4
)
This formula is based on the maximum likelihood estimation for Logistic regression, with an added regularisation term a b, b to ensure that a minimum always exists and avoid overfitting. In this work we always set a = 0.1. The above minimisation problem can be solved by a gradient descent method. Denote by b the solution of the optimisation problem above.

For a new object x new , the probabilistic prediction based on logistic regression will be:

p new = 1 1 + e -x new , b , (5)
which estimates the maximum likelihood probability that y new = 1 if the data are generated by a distribution from a logistic model. We will call p new a direct prediction to distinguish it from multi-probabilistic predictions produced by a Venn machine.

Logistic Regression as an Underlying Algorithm

Now we can describe how logistic regression can be plugged in Venn machine as an underlying algorithm. As earlier the aim is to predict labels y i which are equal to 0 for the controls and 1 for the cases, by objects x i -vectors of features, which are intensities of the most frequent peaks in the logarithm scale.

The probabilistic method of logistic regression allows us to create a new type of taxonomy. The logistic taxonomy τ i , i = 1, . . . , n + 1 is defined as follows. The solution of the optimisation problem b is calculated for the whole set (x 1 , y 1 ), . . . , (x n , y n ), (x n+1 , y) as a training set and is used to make direct predictions p 1 , . . . , p n+1 on the same (training) examples. These predictions are not fair leave-one-out predictions, but it is correct to use them for taxonomy construction.

Let p i , i = 1, . . . , n + 1 be direct predictions p i sorted in the ascending order. Set a number of taxonomy categories K. Each of them will contain approximately equal number of examples. We use K = 5 in this research (the dependence on the parameter K is discussed in Appendix). Let L 0 = 0 and L 1 , L 2 , . . . , L K be the integers closest to (n + 1)/K, 2(n + 1)/K, 3(n + 1)/K, . . . , n + 1. The category τ i is then defined accordingly to the intervals formed by division points p L0 , p L1 , . . . , p L K where value p i falls: τ i is the number j ∈ {1, . . . , K} such that p Lj-1 < p i ≤ p Lj .

Mass Spectrometry Data

We would like to develop a methodology of providing multiprobability predictions for proteomic mass spectrometry data. Hence, we have to take into account the format of the data and peculiarities of the problem.

Mass spectra plots can be noisy because of physical, electrical or chemical sources. Pre-processing is applied to mass spectra to get rid of these systematic artefacts. Auxiliary goals of pre-processing are to normalize the spectra from different samples and reduce the dimensionality of the data. Pre-processing can include the following steps: smoothing by averaging the intensities within a moving window; baseline subtraction; normalization to make sure that the total amounts of ions across different samples are the same. After the true signal is extracted from mass spectra, peaks are identified in each spectrum and then aligned, that is, peaks from different spectra get related to each other and are considered as one peak. Finally, the intensities of identified peaks are calculated. Detailed description of the preprocessing of mass spectrometry samples can be found in [START_REF] Devetyarov | Analysis of serial UKCTOCS-OC data: discriminating abilities of proteomics peaks[END_REF].

Thus, the data we apply in our methodology is represented as intensities of identified peaks. The peaks are usually sorted by their frequency: the more examples a peak is presented, the higher the rank of the peak. We usually consider a certain number of the most frequent peaks only. Thus, every object x i is a vector of features, which are intensities of the most frequent peaks. Each sample of the OC data set is also assigned a level of OC biomarker CA125 that helps discriminate OC samples from healthy samples. For this reason, for OC we will make our diagnosis based not only on MALDI-TOF data, but also on CA125 levels.

The total number of samples is: 561 in HD data set (187 cases and 374 controls); 312 in OC data set (104 cases and 208 controls). Originally, each case was accompanied by two controls matched on patient age, sample collection location and sample collection date/time, among other factors. For this reason, in each data set the number of controls is twice as greater than the number of cases.

Each object, x i , only comprises intensities of the most common peaks. It was shown in statistical analysis [START_REF] Tiss | Highly accurate detection of ovarian cancer using CA125 but limited improvement with serum MALDI-TOF MS profiling[END_REF] that the information useful for discrimination between healthy and diseased samples is concentrated in peaks 2 and 3 in OC data. For this reason, each OC object comprises the five most common peaks. As for the HD data, we consider the peaks represented in at least 1/3 of samples (41 peaks) in this data set.

Each case is assigned a non-negative value T (τ ) -time to diagnosis confirmed by histology/cytology for OC data and time to death for HD data. We will refer to this value as 'time to diagnosis/death'. Each sample in OC data set is also assigned a level C of the biomarker CA125. We will not be trying to predict this time, but we need it to form cross-validation sets for the early diagnostics.

Results and Discussion

To demonstrate how the proposed methodology works in practice, we applied the designed algorithms to HD and OC proteomic data sets.

In all experiments, we use leave-one-out mode: each example (x i , y i ) is considered as if it were a new test example and all the remaining examples in the data are treated as the training set.

We are applying Venn machines with the taxonomy based on logistic regression to MALDI-TOF datasets. Each object x i is a vector comprising the following features: of the most frequent peaks, CA125 value (for OC dataset), additional dimension constantly equal to 1.

Since logistic regression also produces probability distributions, we can compare the results of the application of the Venn machine based on logistic regression and the probabilistic predictor of logistic regression itself. The experiments were applied to the same type of objects x i in the leave-one-out mode. ). As prediction interval indicate, the correct labels for example 1 and 2 are 0 and 1, respectively. The table also includes predictions p new output by logistic regression for each example. Recall that we call these predictions direct predictions as opposed to Venn predictions output by Venn machines. The table demonstrates that both direct and Venn predictions can be correct or erroneous.

First, we would like to demonstrate validity of Venn predictions: true probabilities of label distribution are covered or almost covered by the interval between lower and upper Venn prediction. Since we do not know true probabilities of label distribution, we compare empirical probabilities, that is, mean true labels with the mean direct and Venn predictions. Firstly, the plot demonstrates validity of Venn machine outputs. Secondly, we can see that probability intervals output by Venn machines are narrow (0.025 on average for HD data); hence, they are almost as precise as single probabilities. Finally, Figure 1 demonstrates that probability intervals can be more accurate than single probabilities produced by logistic regression. It can be seen from the Figure that the true labels are very different from the direct predictions but are only slightly above the upper Venn prediction up to approximately 210 examples and within the upper and lower Venn predictions for the remaining examples after this point. Thus, direct predictions can be misleading (the cancer Fig. 1. Cumulative Venn and direct predictions for the HD data probability is oversetimated, while Venn predictions cover true labels: an average true label lies between average low and upper bounds given by the Venn Machine.

It can be said that both algorithms relied on the assumption of the mechanism generating the data -the logistic regression statistical model. However, probability predictions used this mechanism directly, and Venn machines deployed the mechanism when defining the taxonomy. As a result, since the statistical model does not hold true (the opposite can be guaranteed only for artificially generated data), probabilities output by logistic regression are different from empirical probabilities. In contrast, Venn machine's validity was not affected by the fact that the model is not correct. Hence, Venn machine predictions appeared to be more accurate than singleton probability predictions.

Even though the Venn machines and logistic regression produce multiprobability and probability predictions, respectively, their outputs can be interpreted as usual bare predictions (forced predictions). However, we should bare in mind that when we force Venn machines to output a single prediction, they lose their theoretically proven property of validity. In this section we examine the accuracy of forced predictions when we are not able to use advantages of multiprobability predictions.

We can extract forced predictions out of Venn machines and logistic regression the most intuitive way: we classify a new sample as 1 (case) if and only if p new > 0.5 for direct prediction or p + new + p - new > 1 for Venn prediction. This will also allow us to compare accuracy of Venn predictions with the direct predictions.

Given that the aim is to predict the disease as early as possible. For this reason we consider the dynamics of predictive ability of mass spectrometry peaks across the timeline: the accuracy of the proposed methodology on samples in different time slots of the fixed interval (6 months).

Table 2 allows us to compare accuracy of the forced predictions by Venn machine and logistic regression for the OC data. The table demonstrates that Venn machines are comparable with logistic regression in terms of forced prediction accuracy: in time slots close to the moment of diagnosis Venn machine is slightly outperformed by logistic regression, then in months 5-7 they have equal accuracy, and in months 8-11 (time slots we are mostly interested in) Venn machine overperforms logistic regression. Venn machines produce predictions with accuracy higher than 73% up to 10 months in advance of the moment of diagnosis.

For HD we consider the whole dataset rather than dynamics across the timeline, since it is sufficient to predict this disease at any moment to prevent the consequences. The accuracy of the application of Venn machines to the HD data is 69.9%. The accuracy is again comparable with the accuracy of underlying algorithms: 67.9%. 3. Dynamics of prediction intervals output by Venn machines for measurements taken from the same OC case Table 3 shows the dynamics of prediction intervals output by Venn machines for samples 29 and 39. Each row corresponds to a single measurement. Column 2 demonstrates how early in advance this measurement was taken. These samples with multiple measurements illustrate two trends in probability interval change. First, the interval is getting narrower when the moment of diagnosis is approaching, which means that two probability distributions produced by Venn machines are getting closer to each other, and as a result, the overall prediction is getting more precise. This also means that the logistic regression as the underlying model becomes more adequate when the time to diagnostics is closer. Second, the interval is moving towards 1. This implies that we have more trust in our prediction and the prediction is indeed correct.

Conclusion

This paper introduced the methodology of hedging predictions for proteomic mass spectrometry data. We applied the described methodology to the MALDI-TOF data sets and demonstrated how it works. We empirically confirmed the validity of Venn machines and demonstrated that Venn machines can provide narrow probability intervals that are more accurate than the probabilities provided by its underlying algorithm.

Even though Venn machines produce multiprobabilistic predictions, their output can be interpreted as predictions without hedging, similarly to the output of conventional machine learning methods. It was demonstrated that when forced to make single predictions, our methodology provides accuracy similar to the accuracy of the underlying algorithms. As a result, this methodology can provide high accuracy well in advance of the moment of the disease diagnosis.

Figure 1

 1 is a graphical representation of corresponding cumulative results. The horizontal axis shows the number of observed examples. The vertical axis shows the cumulative values of: (1) true labels y new (a solid line); (2) lower and upper Venn predictions p - new , p + new (two dot-dashed lines) and (3) cumulative direct predictions p new (a dashed line). The examples are sorted according to direct predictions.

Table 1 .

 1 Leave-one-out Venn predictions for HD data Results of experiments for several controls and cases of HD data are shown in Table 1 for illustrative purposes. For each example, the table contains the true label y new , and the probability interval [p - new , p + new ]. For example, Venn machines output prediction intervals [0.313, 0.321] and [0.616, 0.616] for probabilities that examples 1 and 2 are cases (y = 1

	No. True label Venn prediction Direct prediction
	1	0	0.313-0.321	0.508
	2	1	0.616-0.616	0.689
	3	0	0.321-0.330	0.510
	4	0	0.143-0.259	0.371
	5	0	0.616-0.634	0.622

Table 2 .

 2 Dynamics of Venn machine and logistic regression performance on the OC dataset

Acknowledgement

This work was supported by EraSysBio+ grant funds from the European Union, BBSRC and BMBF "Living with uninvited guests: comparing plant and animal responses to endocytic invasions" to the Salmonella Host Interactions Project European Consortium; MRC grant G0802594 (Application of conformal predictors to FMRI research); MRC grant G0301107 (Proteomic analysis of the human serum proteome); Veterinary Laboratories Agency (VLA) of Department for Environment, Food and Rural Affairs (Defra) on Machine learning algorithms for analysis of large veterinary datasets; a grant from The National Natural Science Foundation of China (No.61128003); and by grant 'Development of New Venn Prediction Methods for Osteoporosis Risk Assessment' from the Cyprus Research Promotion Foundation.

Appendix. Dependence on the taxonomy parameter

In this work we used K = 5. On Figure 2 we show what happens if this parameter is changed. It can be seen that with the accuracy becomes satisfactory with at least 4-5 taxa. On the other hand when the number K of taxa increases the interval width becomes wider (less informative) without essential improvement of the accuracy. So the choice of K = 5 was reasonable enough.


Stefanos Vrochidis
email: stefanos@iti.gr

Victor Epitropou
email: vepitrop@isag.meng.auth.gr

Anastasios Bassoukos

Sascha Voth
email: sascha.voth@iosb.fraunhofer.de

Kostas Karatzas

Anastasia Moumtzidou
email: moumtzid@iti.gr

Jürgen Mossgraber
email: juergen.mossgraber@iosb.fraunhofer.de

Ioannis Kompatsiaris

Ari Karppinen
email: ari.karppinen@fmi.fi

Jaakko Kukkonen
email: jaakko.kukkonen@fmi.fi

Jürgen Moßgraber

Extraction of Environmental Data from On-Line Environmental Information Sources

Keywords: Environmental, air quality, heatmap, image processing, OCR, data reconstruction, template configuration

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Analysis of environmental information is considered of utmost importance for human population, as this is strongly related to health issues (e.g. cardiovascular diseases), as well as to a variety of important activities (e.g. agriculture). In everyday life, environmental conditions of the atmospheric environment, in terms of air quality, weather, pollen measurements and forecasts are also of particular interest for outdoor activities (e.g. trip planning) and therefore they strongly affect the quality of life. Nowadays, the main sources of such information for the everyday user are web portals and sites. In order to support people in everyday action planning considering the environmental conditions, we need to provide them with services, which combine complementary environmental information from several resources, with a view to generate more reliable environmental measurements. The first step towards this direction is the extraction of data from environmental resources. In practice only a few of the data providers make available some means of access to their actual (numerical) forecast data. In this context, this paper addresses the semi-automatic extraction of air quality forecasts from heatmap images.

After studying a number of on-line chemical weather forecasts by various providers [START_REF] Balk | A European open access chemical weather forecasting portal[END_REF], it can be said that the air quality information is most usually presented in the form of images representing forecast pollutant concentrations over a geographically bounded region, typically in terms of maximum or average air pollution concentration values for the time scale of reference, which is usually the hour or day [START_REF] Karatzas | Internet-based management of Environmental simulation tasks[END_REF], [START_REF] San José | Computational Air Quality Modelling[END_REF], [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF], [START_REF] Karatzas | A European Chemical Weather forecasting Portal. 31st ITM -NATO[END_REF], [6]. These providers present their air quality forecasts almost exclusively in the form of preprocessed images with a color index scale indicating the concentration of pollutants. In addition, these providers arbitrarily choose the resolution of their images, the color scale and color depth employed for visualizing pollution loadings, the covered region, as well as the geographical map projection. The actual mode of presentation varies from simple web images to more elaborated AJAX, Java or Adobe Flash viewers [START_REF] Kukkonen | COST ES0602: Towards a European network on chemical weather forecasting and information systems[END_REF]. While this representation is informative for the casual user (e.g. compared to a table with numerical values), it has the drawback that the data are being presented in a wide range of highly heterogeneous forms, which makes it very complicated to extract and compare their results. To make it worse, some of the images are permanently marked with visible watermarks, text, lines etc. that would make the extraction phase even more challenging.

In order to address this challenge we propose a semi-automatic framework for extracting air quality information from such images and store them into a numerical format. The proposed system is based on an annotation tool, which supports an administrative user to generate a configuration template for each heatmap, and on Optical Character Recognition (OCR) techniques for text information extraction. The basic functionality of the system (i.e. the information extraction from heatmaps), is based on AirMerge [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF], [6], [START_REF] Epitropou | Evaluation of the accuracy of an inverse image-based reconstruction method for chemical weather data[END_REF], [START_REF] Epitropou | Fusion Of Environmental Information For The Delivery Of Orchestrated Services For The Atmospheric Environment In The PESCaDO Project[END_REF], a system that allows for the automatic harvesting, annotation, harmonization and reverse engineering of heatmaps, in order to come up with easily deployable numerical values of chemical weather forecasts.

The contribution of this paper is the methodology and the framework for userassisted air quality information extraction from heatmaps, which extends previous works (i.e. AirMerge) by further adding OCR techniques, as well as allowing user configuration with the aid of a dedicated graphical user interface. More specifically, we propose a framework, which is based on a novel heatmap Annotation Tool (AnT), on the application and optimization of OCR techniques for textual information extraction from heatmap images and on the AirMerge tool [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF] for image processing. This paper is structured as follows: section 2 presents the related work, while section 3 describes the framework architecture. Section 4 presents the Annotation Tool, section 5 the OCR techniques and section 6 the AirMerge system. The results are presented in section 7 and finally, section 8 concludes the paper.

Related Work

Existing maps can be grouped into map types based on the placement and presentation of their information. Discriminating factors between map types can be found in their scale, colorization, quality, accuracy, topology and many other aspects. In case of air quality (or chemical weather) maps there are mainly two types of information covered by the map data: a) Geographical information: points and lines describing country frontiers or other well-known points of interests or structures (e.g. sea, land, capitals) in a given coordinate system, b) Color information: measured data of any kind (e.g. average temperature), which are coded via a color scale representing the measured values. Single values are referenced geographically by a color value at the corresponding geographical point. Chemical weather maps often use this type of maps called raster map or heatmaps images to represent measured data. There are several approaches to extract and digitalize this image information automatically. First, the authors in [START_REF] Musavi | Map processing methods: an automated alternative[END_REF] describe the process of the vectorization of digital image data. Hereby the geographical information, in form of lines, is extracted and converted to digital storable vector data. Only the lines are processed. The work in [11] makes use of the specific knowledge of the known colorization in USGS maps, to have the ability to automatically segment these maps based on their semantic contents (e.g. roads, rivers). In [START_REF] Cao | Text/graphics separation in maps[END_REF] the segmentation quality of text and graphics in color map images is improved, to enhance the results of the following analysis processes (e.g. OCR), by selecting black or dark pixels from color maps, cleaning them up from possible errors or known unwanted structures (e.g. dashed lines), to get cleaner text structures.

Although research work has been conducted towards the automatic extraction of information in maps, very few works address the automatic extraction of information from chemical weather maps. In such works [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF], [6], [START_REF] Epitropou | Evaluation of the accuracy of an inverse image-based reconstruction method for chemical weather data[END_REF], a method to reconstruct environmental data out of chemical weather images is proposed. In a first step the relevant map section is scraped from the chemical weather image. After that disturbances are removed (e.g. country lines) and a color classification is employed to classify every single data point (pixel), to recover the measured data. With the aid of the known geographical boundaries, given by the coordinate axis and the map projection type, the geographical position of the measured data point can be retrieved. In case of missing data points, a special interpolation algorithm is used to fill these gaps.

The proposed work goes one step beyond the aforementioned heatmap extraction methods, since it introduces a configurable user-assisted environment, which facilitates the application of the framework on new heatmaps without requiring programming skills and low level configuration on the user's side.

Framework Architecture

The architecture of the proposed framework is illustrated in Figure 1 and includes two main components: the Annotation Tool and the data extraction service.

The first phase, called "Template Configuration" (1→2→3), includes the manual annotation of an image with the AnT, and the generation of a configuration file. This process is controlled by an administrative user with the aid of AnT. The second phase includes the "Data extraction" (1+3→4→5), which uses the configuration file to extract data from the specific heatmap. During this phase, the parts of each image are analyzed using image and text processing techniques. Specifically, the heatmap is processed with the AirMerge system, while the text information located in the image is extracted and processed using OCR techniques and text processing.

Annotation Tool

The Annotation Tool (AnT) is used to interactively annotate heat maps and it was developed in order to make the annotation process easier for the user.

To make the tool platform independent, the QT Framework1 was used. The implementation is designed via the MVC (Model/View/Controller) pattern, to ensure its expandability. To allow for different interaction possibilities, two views were implemented. First a simple Tree View, which represents the XML structure and its entries as traversable tree, and a window, which represents the selected tree data graphically. Regions of Interest (ROIs) and Points of Interest (POIs) are drawn onto this window.

OCR techniques

The OCR module uses the information of the configuration file to extract textual data from images and improve the results using text processing based on heuristic rules.

The first processing step of this module includes the application of OCR on specific parts of the input image as the title, the color scale, the map x and y axes parts. The OCR software that was used is Abbyy Fine Reader 3 .

In the second step, we apply text processing based on heuristic rules in order to correct, extract and understand the semantic information encoded in the aforementioned locations. Each of these locations was treated in a different way.

The title, if it exists, usually contains the name of the aspect, the measurement units and the date/time. The measurement units are usually standard depending on the measured aspect so we do not need to extract them. The date/time is considered as the element that is the most difficult to extract, given the fact that many different formats exist. In order to correct possible mistakes in the textual format of the month, day or aspect we exploited the Levenshtein distance. More specifically, three English ground truth sets were created for the three aforementioned elements and were compared to the corresponding OCR result. Then, we have selected the element from the ground truth that had the smaller Levenshtein distance from the text generated by the OCR.

The color scale contains the values that each color of the map corresponds to. The processing and extraction of information from the color scale element can be divided into two main parts. In the first, we attempt to check and correct OCR results for the scale, while in the second we correlate values to colors. In order to correct the OCR results, we find the most common range among the scale values and adapt accordingly the mistaken values. The correlation of values to colors is achieved by pointing at the middle of each color by using the coordinates of the values.

The last two elements that are analyzed with OCR are the x and y axes. In the case of heatmaps the two axes contain similar information and thus we will apply similar processing techniques to them. The information that can be obtained by each axis is the geographical coordinates of the points of the map. In order to realize this, we have to identify successfully at least two points (x, y coordinates) of the map axes and define their position in relation to the map. In order to identify these points, we first correct most of the errors produced by OCR and then use the coordinates of the elements, as defined by OCR to specify the position of those points.

AirMerge

AirMerge is a web-based system that supports harvesting Chemical Weather forecast images and converting them to numerical data. A derivative of its image processing engine is used in the Data Extraction phase of the proposed toolchain, and it is already used for creating harmonized, numerical Chemical Weather data 4 .

The AirMerge system combines elements of screen scraping and innovative image processing algorithms [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF], [6], [START_REF] Epitropou | Evaluation of the accuracy of an inverse image-based reconstruction method for chemical weather data[END_REF] in order to produce uniform, indexed data. These data are then stored in a back-end database and may be recalled for further processing such as numerical applications, model ensembles, visualization, transformation etc. The main task of AirMerge is the extraction of data from heatmaps. This is achieved by using a processing chain that consists of a "screen scraping" phase, where raw RGB pixel data are extracted from heatmaps, a mapping phase, where RGB values are classified to a color scale and mapped to ranges of numerical values and a linear deprojection phase, where the images' raster is interpreted as a geographical grid in a specified geographical projection, centered on key points. Screen scraping procedure: This step handles the cropping of the original image to a region of interest and parsing of it into a 2D data array directly mapped to the original images' pixels. Also, it associates the color to minimum/maximum value ranges of the air pollutant concentration levels, which is often implied by the color scale associated with the original images. It should be noted that the information about where to crop, where each color on the legend is, to which index it should correspond, etc. are provided by the configuration template of the AnT in the proposed system. In this phase, the mapping of the images' raster to a specific geographical grid is per-formed, since the images themselves represent geographical region. The configuration system allows choosing between the most commonly encountered geographical projections (equirectangular, conical, polar stereographic etc.) and choosing keypoint in the image to allow for precise pixel-coordinate mapping.

"Reconstruction of missing values and data gaps" procedure: This step is introduced to deal with unwanted elements such as legends, text, geomarkings and watermarks, as well as regions that are not part of the forecast area, which might be present after the screen scraping phase. The image pixels are classified into three main categories: valid data (with colors that satisfy the color scale's classification), invalid data (with colors not present in the color scale), and regions containing colors that are explicitly marked for exclusion, and which are considered void during processing. Such marked regions are not considered as part of the forecast, and thus do not undergo data correction. However, regions containing unmarked invalid data are considered as regions with correctable errors or "data gaps" which can be filled-in. This distinction is due to their different appearance patterns: void regions are usually extended and continuous (e.g. sea regions not covered by the forecast, but present on the map), while invalid data regions are usually smaller but more noticeable (e.g. lines, text, watermarks etc.) and with more noise-like patterns, and thus it is more compelling to remove them by using gap-filling techniques. These techniques include traditional grid and pattern-based interpolation techniques using neural networks.

It should be noted that the AirMerge system functionality is also provided via an API, which is available as a REST service [START_REF] Epitropou | Fusion Of Environmental Information For The Delivery Of Orchestrated Services For The Atmospheric Environment In The PESCaDO Project[END_REF]. Therefore, AirMerge can serve any request related to the heatmaps of many chemical weather models (e.g. every-day processing), thus making it suitable for environmental service-oriented applications.

Results

In this section, we present the results of the framework when applied in three air quality heatmaps from different providers. Since an evaluation of AirMerge is already provided in [START_REF] Epitropou | Evaluation of the accuracy of an inverse image-based reconstruction method for chemical weather data[END_REF], we evaluate the results of the OCR and the total system output. Regarding the OCR, we focus on the recognition of the x and y axes, since this is the most important information in order to correctly map the air quality index onto the right coordinates. The following providers are considered for the evaluation: GEMS 5 , Laboratory of Atmospheric Physics of the Aristotle University of Thessaloniki 6 and the Atmospheric and Oceanic Physics Group7 .

GEMs website

Figures 3 and4 depict the original and reconstructed image by the Airmerge system, which are almost identical, and any noise (e.g. black lines) was removed [START_REF] Epitropou | A method for the inverse reconstruction of environmental data applicable at the Chemical Weather portal[END_REF], [START_REF] Epitropou | Evaluation of the accuracy of an inverse image-based reconstruction method for chemical weather data[END_REF]. During this process an error is usually introduced mostly due to the inability of OCR to perfectly identify the position of each coordinate on the map axes. In table 1 we report the longitude and latitude steps (i.e. the coordinate step for each pixel), the correct step value between two subsequent coordinate marks (e.g. when the marks are 0 0 , 5 0 , 10 0 , etc. the step value is 5), the estimated value, the absolute difference and finally the introduced error. In both cases the error is very low and acceptable (in general we assume that an error is acceptable, when it is less than 3%).

Laboratory of Atmospheric Physics of the AUTH site

In a similar way we present the initial and the reconstructed image of this website in Figures 5 and6. The results are reported in table 2 and the error is again very small.

Atmospheric and Oceanic Physics Group site

Finally, in table 3 we present results for the last provider reporting an average error of 0.35%. The initial and the reconstructed map are illustrated in figures 7 and 8. It should be noted that the white region in figure 7 is treated as "void space" in figure 8, and considered as a distinct case than national border lines, which are instead treated as unwanted noise and filled-in.

Conclusions

Despite the fact that the current presentation format of air quality forecasts might be ideal for casual users, it is not easily accessible by automatic services which would expect a structured and numerical format of the forecast data. In this context, we have proposed a framework for air quality information extraction from heatmaps, combining existing (AirMerge), as well as new (AnT and OCR) components. This framework could serve as a basis for supporting environmental systems that provide either air quality information from several providers for comparison or orchestration purposes or high level suggestions on everyday issues (e.g. travel planning) based on advanced decision support [START_REF] Wanner | Building an Environmental Information System for Personalized Content Delivery[END_REF], which could facilitate the quality of life. The proposed work overcomes the limitation of not having access to the raw data, since it only considers information being publicly available on the Internet, thus respecting also data access policies. Although the system has been tested with forecast air quality heatmaps it could also deal with observed pollutant and pollen concentrations represented in the same way. Future work includes extensive evaluation with more images in different projections (e.g. conical) and addressing of pollen heatmaps.

Fig. 1 .

 1 Fig. 1. Air quality data extraction framework.

Figure 2

 2 depicts the AnT tool after a heatmap from GEM's Project 2 site is loaded. The air quality heatmaps contained in the site are typical examples of images used for representing chemical weather forecasts. The left part of the tool contains the heatmap as well as the ROIs, which are the following: a) the map itself, b) the x and y axis related to the heatmap, c) the color scale, d) the numbers corresponding to the color scale and e) the title of the heatmap. The ROIs are depicted as red bounding boxes and are defined by the user. Finally, their values are recorded to the right part of the AnT inside the XML template under the corresponding nodes.

Fig. 2 .

 2 Fig. 2. Screenshot of Image Annotation tool

Fig. 3 .Fig. 4 .

 34 Fig. 3. Original Image

Fig. 7 .Fig. 8 .

 78 Fig. 7. Original Image

Table 1 .

 1 Results for the GEMs website ResolutionSteps

			Correct Value	Estimated Value	Absolute Dif-ference	Error
	Longitude step	0.0791	5	5.0634	0.0634	1.25%
	Latitude step	0.0777	5	4.9776	0.0224	0.45%

Table 2 .

 2 Results for the Atmospheric and Oceanic Physics Group website

		Resolution Steps	Correct Value	Estimated Value	Absolute Dif-ference	Error
	Longitude step	0.0311	2	1.9924	0.0076	0.4%
	Latitude step	0.0276	1	1.0236	0.0236	2.3%

Table 3 .

 3 Results for the Atmospheric and Oceanic Physics Group website ResolutionSteps

			Correct Value	Estimated Value	Absolute Dif-ference	Error
	Longitude step	0.0289	2	1.9937	0.0063	0.3%
	Latitude step	0.0249	1	0.9958	0.0042	0.4%

http://qt.nokia.com/products/

http://gems.ecmwf.int/d/products/raq/

http://www.abbyy.com/ Heatmap and component parts XML Template

http://projects.isag.meng.auth.gr/airmerge/

http://gems.ecmwf.int/d/products/raq/

http://lap.physics.auth.gr/forecasting/fore_images/.

http://www.fisica.unige.it/atmosfera/bolchem/MAPS/

Acknowledgments. This work was supported by PESCaDO project (FP7-248594).


A Panayotopoulos

P Vlamos
email: vlamos@ionio.gr

Cutting degree of meanders

Keywords: Dyck path, Grand Dyck path, 2-Motzkin word, elevated peakless Motzkin path

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The cutting problem of meanders has been introduced by the authors in [START_REF] Panayotopoulos | Meandric polygons[END_REF]. Using meandric polygons, we have also enumerated the set of meanders with given numbers of cuttings, i.e., the numbers of arcs of the meandric curve lying above and below the horizontal line at a given point.

In the present paper, we complete the cutting problem of meanders. At first, we represent a meander or meandric system using Motzkin words, and we specify the subset of meandric Motzkin words, through a necessary and sufficient condition. Additionally, we determine the matching nested sets of a meandric Motzkin word. Afterwards, we establish a bijection between meandric Motzkin words and elevated peakless Motzkin paths. We calculate the cardinality of this class, as well as the number of these paths, passing through a given point. Finally, we confront the two main cutting problems of meanders, as we determine the set of meanders with a given sequence of cutting degrees, or with a given cutting degree.

The two cutting problems can be considered as data-mining problems, since we extract the subset of Motzkin words which are meandric, or the subset of Motzkin paths which are peakless and elevated. Specifically in the field of Bioinformatics, the already mentioned problems have many facets in relation with the M-folding problem and RNA prediction problem, based on the fact that closed RNA structures can be viewed as mathematical objects obtained by abstracting topologically non-relevant properties of planar folding of single-stranded nucleic acids. In the literature, many approximation data-mining algorithms had been constructed, concerning the enumeration of sequential terms of RNA secondary structures, through the bijection with alternative representations such as energy models, plane trees and Motzkin numbers [START_REF] Stein | On some new sequences generalizing the Catalan and Motzkin numbers[END_REF], k-noncrossing set of partitions [START_REF] Jiang | A general edit distance between RNA structures[END_REF], set of nested arcs [START_REF] Jiang | A general edit distance between RNA structures[END_REF], Motzkin paths and Dyck paths [START_REF] Nebel | Combinatorial Properties of RNA Secondary Structures[END_REF].

The following definitions and notation refer to notions that are necessary for the development of the paper. A set S of disjoint pairs of [2n] such that {a,b}∈S {a, b} = [2n] and for any {a, b}, {c, d} we never have a < c < b < d, is called nested set of pairs on [2n]. We denote the set of all nested sets of pairs on

[2n] by N 2n . Two nested sets S 1 , S 2 ∈ N 2n define a permutation σ on [2n], such that σ(2i -1) = j iff {2i -1, j} ∈ S 1 and σ(2i) = j iff {2i, j} ∈ S 2 , for every i ∈ [n].

Meanders

A (closed) meander of order n is a closed self avoiding curve, crossing a horizontal line 2n times [START_REF] Lando | Plane and projective meanders[END_REF]. We denote by M 2n the set of meanders of order n [START_REF] Panayotopoulos | Properties of meanders[END_REF].

We can extend the definition of meanders to meandric systems with k components, by allowing configurations with k disconnected meanders [START_REF] Lando | Plane and projective meanders[END_REF].

In this paper, every meander of M 2n is represented by a word w ∈ {a, ā, b, b} * , which is called meandric Motzkin word [START_REF] Panayotopoulos | Meanders and Motzkin words[END_REF]. For example, the meander of In order to formally introduce the meandric Motzkin words, following [START_REF] Panayotopoulos | Meanders and Motzkin words[END_REF], we recall that for every pair

v = v 1 v 2 • • • v ℓ and v ′ = v ′ 1 v ′ 2 • • • v ′ ℓ of words in {a, ā} * of length ℓ, we define v • v ′ to be the word w = w 1 w 2 • • • w ℓ ∈ {a, ā, b, b} * , with w i = v i • v ′ i , i ∈ [ℓ], where a • a = a, ā • ā = ā, a • ā = b, ā • a = b.
Note that the words v and v ′ can be uniquely recovered from w using the following relations

v i = a, if w i = a or b, ā, if w i = ā or b, and v ′ i = a, if w i = a or b, ā, if w i = ā or b, for i ∈ [ℓ]
, and we call v, v ′ associates of w.

We define the set of the meandric Motzkin words of length 2n as the set

Ŵ2n = {w : w = v • v ′ , v, v ′ ∈ D 2n }.
Notice that from the word v • v ′ we immediately obtain the word v ′ • v, by interchanging the letters b and b. So, in order to generate the set Ŵ2n it is actually enough to construct half of its elements.

The set Ŵ2n can be partitioned into the classes

W k 2n , k ∈ [n]
, where W 2n = W 1 2n . Each class W k 2n contains the meandric Motzkin words with v, v ′ ∈ D 2n , such that their corresponding nested sets are k-matching.

It has been proved in [START_REF] Panayotopoulos | Meanders and Motzkin words[END_REF] that there exists a bijection between the sets M 2n and W 2n . According to this bijection, every meander µ can be considered as an ordered pair (U, L) of nested sets of pairs, where U (resp. L) is defined by the arcs of the meandric curve above (resp. below) the horizontal line. Since U, L are encoded by two Dyck words, say v and v ′ , then the meandric Motzkin word w = v • v ′ encodes the meander µ.

In [START_REF] Panayotopoulos | Meanders and Motzkin words[END_REF], it was also proved (see Proposition Proof. Let w ∈ {a, ā, b, b} * and let w = pq be a factorization of w. Then, the associates v, v ′ of w, such that w = v • v ′ , are uniquely defined and, furthermore, the factorization of w uniquely defines the factorizations v = rs and v ′ = r ′ s ′ , such that p = r • r ′ and q = s • s ′ . By definition, we have the following relations

|r| a = |p| a + |p| b , |r| ā = |p| ā + |p|b, |r ′ | a = |p| a + |p|b, |r ′ | ā = |p| ā + |p| b , which lead to |r| a -|r| ā = |p| a -|p| ā+|p| b -|p|b and |r ′ | a -|r ′ | ā = |p| a -|p| ā-(|p| b -|p|b). (1)
Note that the above relations also hold, when we set p = w, and in this case we have r = v and r

′ = v ′ . If w ∈ Ŵ2n , then v, v ′ ∈ D 2n , so that |r| a -|r| ā ≥ 0, |r ′ | a -|r ′ | ā ≥ 0 and |v| a -|v| ā = |v ′ | a -|v ′ | ā = 0.
-|v| ā = |v ′ | a -|v ′ | ā = 0. From the above, we deduce that v, v ′ ∈ D 2n , so that w ∈ Ŵ2n .
The process of determining the matching nested sets U, L that correspond to a meander w ∈ W 2n can be simplified by the following easily verified proposition: Proposition 2 Let U, L be the matching nested sets that correspond to a meander w ∈ W 2n . Then, every pair of U (resp. L) corresponds to a subsequence of w having one of the following forms: aā, a b, bā, b b (resp. aā, ab, bā, bb).

According to Proposition 2, the short pairs of U (i.e., pairs of the form {i, i + 1}) are immediately determined by identifying the positions of the occurrences of the digrams aā, a b, bā, b b in w. Then, by deleting these digrams from w, we obtain a new word w ′ ∈ {a, ā, b, b} * . We repeat the same procedure for w ′ , until all pairs are determined. The pairs of L are obtained analogously.

For example, for the meander w = aa bbabāā bā of Fig. 1, the nested set U is obtained as follows: and by the fact that its corresponding nested sets

U = {{1, 10}, {2, 7}, {3, 6}, {4, 5}, {8, 9}}, L = {{1, 4}, {2, 3}, {5, 8}, {6, 7}, {9, 10}},
which are obtained using Proposition 2, are 2-matching.

Motzkin paths

We recall that a Motzkin path (resp. ). An elevated peakless Motzkin path, is a Motzkin path that never touches the x-axis except at its first and last points and it has no peaks (i.e., no occurrences aā). For example, the path of Fig. 2 peakless Motzkin paths of length 2n.

Given that (see [START_REF] Sloane | The Online-Encyclopedia of Integer Sequences[END_REF], A004148) the number of peakless Motzkin paths of length n is equal to a

(n) = k≥⌈ n+1 2 ⌉ 1 k k n-k k n-k+1
, we easily deduce that the number of paths of P 2n is equal to a(2n -2). Thus, we have the following result.

(i) if a vertex labeled by ξ belongs to the e th level, then ξ ≤ min{e, 2n-e}, (ii) to every vertex with label ξ correspond at most three children, with labels ξ -1, ξ, ξ + 1, on [n]. (iii) at every path of the tree, we can not find three consecutive vertices with labels ξ, ξ + 1, ξ on [n].

Our next effort concerns the number a i,k of paths of P 2n passing through the point (i, k). We already know (see [START_REF] Sloane | The Online-Encyclopedia of Integer Sequences[END_REF], A097724) that the number of left factors of peakless Motzkin paths, having length n and endpoint height k is equal to

T (n, k) = (k + 1) j≥⌈ n-k+1 2 ⌉ 1 j j n-k-j j+k n+1-j .
The number of left factors for the paths of P 2n starting from (0, 0) and ending at (i, k) is given by

T (i -1, k -1) = k j≥⌈ i-k+1 2 ⌉ 1 j j i-k-j j+k-1 i-j
. Furthermore, the number of right factors for the paths of P 2n starting from (i, k) and ending at (2n, 0) is equal to the number of left factors for the paths of P 2n starting from (0, 0) and ending at (2ni, k), i.e., T (2n

-i -1, k -1) = k j≥⌈ 2n-i-k+1 2 ⌉ 1 j j 2n-i-k-j j+k-1
2n-i-j . Thus, we have the following result.

Proposition 5

The number of paths of P 2n passing through the point (i, k) is given by the product k 2

j≥⌈ i-k+1 2 ⌉ 1 j j i-k-j j+k-1 i-j j≥⌈ 2n-i-k+1 2 ⌉ 1 j j 2n-i-k-j j+k-1 2n-i-j .
The paths of P 2n passing through the point (i, k) are obtained by the subtree of the tree of H 2n-1 with η(i) = k.

We denote by P

2n the set of elevated peakless 2-Motzkin paths y of length 2n, with |y| b = |y|b. Since the length of every path w ∈ P 2n and the number |w| a + |w| ā are even, we conclude that |w| b is also even. Therefore, every path y ∈ P Furthermore, the number of peakless Motzkin paths of length 2n with k a steps and 2ν horizontal steps is given by the formula (sequence A089732 in [START_REF] Sloane | The Online-Encyclopedia of Integer Sequences[END_REF])

T (2n, n -ν) = 2n-k k 2n-k k+1 1 2n-k = n+ν 2ν n+ν 2ν-1 1
n+ν , where T (2n, nν) is also the number of paths in P 2n with 2ν horizontal steps, for n, ν > 0.

Finally, since P 2n is partitioned with respect to the number 2ν of b steps, where 1 ≤ ν ≤ n, using the previous formula, we have the following result.

Cutting problem

For every meander of M 2n and for any i ∈ [2n -1] we consider the vertical line passing through the middle point of the segment (i, i + 1), which we call as the i-line [START_REF] Panayotopoulos | Meandric polygons[END_REF]. The number of arcs of the meandric curve, which are intersected by the i-line is called cutting degree ǫ(i) of the meander at i. For example, for the meander of Fig. 1, we have ǫ = 244466422 (see Fig. 3).

The cutting degrees for a meander w ∈ W 2n are obviously deduced by

ǫ(1) = ǫ(2n -1) = 2 and ǫ(i) =      ǫ(i -1) + 2, if w i = a, ǫ(i -1), if w i = b or b, ǫ(i -1) -2, , if w i = ā, for i = 2, 3, . . . , 2n -2.
We denote by E 2n-1 the set of the sequences ǫ = (ǫ(i)

) i∈[2n-1] satisfying ǫ(1) = ǫ(2n -1) = 2 and |ǫ(i) -ǫ(i + 1)| ∈ {0, 2, -2}, i ∈ [2n -2].
From the previous, we have that 1. To every meander w ∈ W 2n corresponds a unique finite sequence ǫ ∈ E 2n-1 but obviously, the converse does not hold. 2. The sequences ǫ and η, corresponding to a meander of M 2n , satisfy ǫ(i) = 2η(i), i ∈ [2n -1]. 3. For the cutting degree, the following properties hold

(i) 2 ≤ ǫ(i) ≤ 2 min{i, 2n -i}, i ∈ [2n -1]. (ii) |ǫ(i + 1) -ǫ(i)| ∈ {0, 2}, i ∈ [2n -1]. (iii) There exists no i ∈ [2n -2], such that ǫ(i) = ǫ(i + 2) = ǫ(i + 1) -2.
The main cutting problems of meanders are to determine:

a) The subset M 2n (ǫ) of M 2n with given sequence of cutting degrees ǫ.

b) The subset M 2n (i, ǫ(i)) of M 2n with given cutting degree ǫ(i) at i.

For the first problem, i.e. to determine the set W 2n (ǫ) we have that if ǫ is deduced by a given meander then W 2n (ǫ) = ∅, while if ǫ is deduced by a given η ∈ H 2n-1 then W 2n (ǫ) could be empty.

So we start by finding the Motzkin word y ∈ P 2n which corresponds to the sequence η ∈ H 2n-1 , where According to the above, the problem can be confronted by the following procedure:

η(i) = ǫ(i) 2 , for i ∈ [2n -
Step 1. Find the Motzkin word y which corresponds to a given sequence η ∈ H 2n-1 .

Step 2. Find the tree B 2ν , where ν = 1 2 |y| b . Step 3. Find, using y and B 2ν , the 2-Motzkin words

ω = ω 1 ω 2 • • • ω 2n of {a, ā, b, b} *
for which:

(i) ω i =      a, if y i = u, b or b, if y i = h, ā, if y i = d, for i ∈ [2n]. (ii) |ω| b = |ω|b = ν.
Step 4. Find the associate words r, r ′ for each ω and check if there are Dyck words.

Step 5. Find from the remaining ω, the ones that satisfy the matching property for the corresponding nested sets.

Examples.

The following examples illustrate the above procedure. We will now solve the second problem, i.e., we will determine the elements of the set W 2n (i, ǫ(i)).

For the existence of such meanders, the conditions that ǫ(i) is even and 2 ≤ ǫ(i) ≤ 2 min{i, 2n -i}, for i = 2, 3, . . . , 2n -2, must hold.

If the cutting degree at i for a meander is ǫ(i), then its corresponding path of P 2n passes through the point (i, η(i)) = i, ǫ(i)

2

. It is natural to find the paths of P 2n passing through the point (i, η(i)) and to check which of them satisfy the matching property.

From Proposition 5 we know the number of paths of P 2n passing through the point (i, η(i)) = i, ǫ(i)

2

. These paths of P 2n can be determined by the subtree of the tree of H 2n-1 presented in Section 2.

For example, for n = 3 and ǫ(4) = 2 the number of paths of P 6 passing through the point (4, 1) is equal to 2. Using the corresponding subtree, we obtain the sequences 1 1 1 1 1 and 1 2 2 1 1 hence, Using the above algorithms, we have obtained the following enumeration results: A192927, A207851, A203019, A208357, A208358 in [START_REF] Sloane | The Online-Encyclopedia of Integer Sequences[END_REF].

  Fig. 1 corresponds to the word w = a a b b a b ā ā b ā.

Fig. 1 .

 1 Fig. 1. A meander with w = a a b b a b ā ā b ā.

2 . 1 )Proposition 1

 211 that every word w of Ŵ2n is a 2-Motzkin word of {a, ā, b, b} * , with |w| b = |w|b. More specifically, if w ∈ Ŵ2n , then by deleting the letters b, b (resp. a, ā) in w, a Dyck word (resp. Grand Dyck word) is obtained, or equivalently, we have that |w| a = |w| ā, |w| b = |w|b and |p| a ≥ |p| ā, for every factorization w = pq. Consequently, every word w ∈ Ŵ2n can be considered as the shuffle of a Dyck word in {a, ā} * and a Grand Dyck word in {b, b} * , where the first and last letter of w are a and ā respectively. Obviously, not every word w satisfying the above conditions encodes a meander, or a meandric system. As an example, consider the word a b b b b ā. The next result states the necessary and sufficient condition for a word to be a meandric Motzkin word. For a word w ∈ {a, ā, b, b} * of length 2n, we have that w ∈ Ŵ2n if and only if |w| a -|w| ā = |w| b -|w|b = 0 and |p| a -|p| ā ≥ |p| b -|p|b , for every factorization w = pq.

  In this case, from relations (1), we have that |p| a -|p| ā ≥ |p|b-|p| b and |p| a -|p| ā ≥ |p| b -|p|b, or equivalently |p| a -|p| ā ≥ |p| b -|p|b , while applying relations (1), for p = w, we obtain that |w| a -|w| ā = |w|b -|w| b and |w| a -|w| ā = |w| b -|w|b, or equivalently |w| a -|w| ā = |w| b -|w|b = 0. Conversely, if |w| a -|w| ā = |w| b -|w|b = 0 and |p| a -|p| ā ≥ |p| b -|p|b , then relations (1) imply that |r| a -|r| ā ≥ 0 and |r ′ | a -|r ′ | ā ≥ 0, while applying relations (1) for p = w, we obtain that |v| a

a 1 a 2 b 3 b 4 a 5 b 6 ā 7 ā 8 b 9 ā 10 -→ a 1 b 4 a 5 ā 8 b 9 ā 10 -→ a 1 b 4 b 9 ā 10 -→ a 1 ā 10 Thus,Examples 1 . 1 |p| b -|p|b 0 1 2 2.b b a b ā b b b ā |p| a -|p| ā 1

 145101410110101121 U = {{1, 10}, {2, 3}, {4, 9}, {5, 8}, {6, 7}}. Analogously, we obtain that L = {{1, 8}, {2, 7}, {3, 4}, {5, 6}, {9, 10}}. The meandric Motzkin word w = a b b b b a b ā b ā, does not belong to the set Ŵ10 . This is verified by the following table: a b b b b a b ā b ā |p| a -|p| ā 1 1 The meandric Motzkin word w = a b b a b ā b b b ā, belongs to the set w ∈ W 10 . This is verified by the following table a fact that its corresponding nested sets U = {{1, 10}, {2, 3}, {4, 7}, {5, 6}, {8, 9}}, L = {{1, 2}, {3, 6}, {4, 5}, {7, 8}, {9, 10}}, which are obtained using Proposition 2, are matching. 3. The meandric Motzkin word w = a a b b b b ā b b ā, belongs to the set w ∈ W 2 10 . This is verified by the following table a a b b b b ā b b ā |p| a -|p| ā 1

2 -

 2 Motzkin path) is a lattice path in the first quadrant beginning at the origin, ending at the x-axis and consisting of the steps a = (1, 1) as the up-step, ā = (1, -1) as the down-step and b = (1, 0) as the horizontal step (resp. a, ā and two kinds of horizontal steps denoted by b and b). Every Motzkin path (resp. 2-Motzkin path) is encoded by its sequence of steps, i.e., a word y ∈ {a, ā, b} * (resp. y ∈ {a, ā, b, b} * ). The set of all Motzkin paths (resp. 2-Motzkin paths) of length n is enumerated by the n-th Motzkin number M n = k≥0 n 2k C k[START_REF] Donaghey | Motzkin numbers[END_REF], (resp. the number C n+1[START_REF] Delest | Algebraic languages and polyominoes enumeration[END_REF]

Fig. 2 .

 2 Fig. 2. The elevated peakless Motzkin path with y = a a b b a b ā ā b ā.

  |y| b = |y|b = ν, is uniquely obtained from a path w ∈ P 2n , with |w| b = 2ν, by converting half of the b steps of w into b steps. Conversely, w ∈ P 2n generates 2ν ν paths of P (2) 2n , which cannot be generated by any other path of P 2n .

Proposition 6

 6 The cardinality of P

Fig. 3 .

 3 Fig. 3. A meander with ǫ = 2 4 4 4 6 6 4 2 2.

  1]. From y we determine all the 2-Motzkin words ω ∈ {a, ā, b, b} * with the following properties: The letter h is replaced by the letters b, b such that |ω| b = |ω|b. Obviously, if |y| b = 2ν there exist 2ν ν such words ω. Those words can be obtained, if we properly place the b's and b's in the positions of the b's, with the help of a binary tree B 2ν of height 2ν. In this tree, each left child is labeled by b and each right child is labeled by b, given that the equal number of b's and b's is ν. Thus, the number of possible meanders belongs to the set {0, 2, 4, . . . , 2ν ν }. We notice that the paths of the right subtree of B 2ν can be produced by replacing the b's with b's at the paths of the left subtree. For example, for the Motzkin word y = a a b b a b ā ā b ā of Fig. 2, we have ν = 2 and the following 4 2 = 6 words: a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, which are obtained with the aid of the tree B 4 .

1 ) 2 6 3 =a b b a b ā b b b ā, a b b a b ā b b b ā}. 3 ) 4 2

 12334 For η = 1 1 1 1 1 2 2 1 1 we have y = a b b b b a b ā b ā and ν = 3. From the left subtree of B 6 , we derive 1 10 words. We can easily observe that none of these words satisfies the condition of the fourth step. Hence, for ǫ = 2 2 2 2 2 4 4 2 2, we have that W 10 (ǫ) = ∅. 2) For η = 1 1 1 2 2 1 1 1 1 we have that y = a b b a b ā b b b ā and ν = 3, hence from the left subtree of B 6 , we derive 10 words, only four of them satisfying the conditions of the fourth and fifth step. Hence, for ǫ = 2 2 2 4 4 2 2 2 2, we have that W 10 (ǫ) = {a b b a b ā b b b ā, a b b a b ā b b b ā, a b b a b ā b b b ā, a b b a b ā b b b ā, a b b a b ā b b b ā, a b b a b ā b b b ā, For the meander of Fig. 1 we have η = 1 2 2 2 3 3 2 1 1, y = a a b b a b ā ā b ā, and ν = 2, hence from the left subtree of B 4 , we derive 1 2 = 3 words. We can easily observe that all these words satisfy the conditions of the fourth and fifth steps. Hence, for ǫ = 2 4 4 4 6 6 4 2 2, we have that W 10 (ǫ) = {a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā, a a b b a b ā ā b ā}.

W 6

 6 (4; 2) = W 6 (2 2 2 2 2) ∪ W 6 (2 4 4 2 2) = = {a b b b b ā, a b b b b ā, a b b b b ā, a b b b b ā} ∪ {a a b ā b ā, a a b ā b ā} = {a b b b b ā, a b b b b ā, a b b b b ā, a b b b b ā, a a b ā b ā, a a b ā b ā}

  The sets S 1 , S 2 are k-matching if and only if σ has k cycles. In the case where k = 1, S 1 , S 2 are simply called matching. A word v ∈ {a, ā} * is called Dyck word if |v| a = |v| ā and for every factorization v = pq we have |p| a ≥ |p| ā where |v| a , |p| a (resp. |v| ā, |p| ā) denote the number of occurrences of a (resp. ā) in the words v and p. A word v ∈ {a, ā} * , such that |v| a = |v| ā is called Grand Dyck word. We denote by D 2n the set of all Dyck words of length 2n and it is well known that the cardinality of D 2n is equal to the n-th Catalan number C n = 1

	n+1	2n n .

Proposition 3 The cardinality of P 2n is given by

A list of terms of the above sequence are presented in [START_REF] Sloane | The Online-Encyclopedia of Integer Sequences[END_REF], A203019.

The set P 2n can also be encoded in {0, 1, 2, . . . , n} by the heights η(i) of the vertices (i, η(i)), i = 0, 1, . . . , 2n of the path. For the heights η(i) we have the following properties:

Indeed, the first property holds since each path of P 2n begins at (0, 0), ends at (2n, 0) and it is elevated. The second property is deduced by the fact that the difference of the number of occurrences of a's and ā's at the i-step is at most i, if i ≤ n or 2ni, if i ≥ n. The third property is obvious since at each step we add +1, -1 or 0 to the height of the previous point. The last property is an immediate consequence of the fact that the path is peakless.

By omitting the initial and final terms η(0) = η(2n) = 0, we introduce the finite sequences η = (η(i)) i∈[2n-1] , and their set will be denoted by H 2n-1 .

For example, for the path of Fig. 2, we obtain the sequence

In general, to every path y ∈ P 2n corresponds a sequence η ∈ H 2n-1 with

The converse also holds, since to every sequence η ∈ H 2n-1 corresponds a path y ∈ P 2n with y 1 = a, y 2n = ā and

So, we have arrived at the following proposition.

Proposition 4 There exists a bijection between the sets P 2n and H 2n-1 .

From the above bijection between the sets P 2n and H 2n-1 , we deduce that the cardinality of the set H 2n-1 is given by the formula of Proposition 3.

The paths of P 2n can be obtained by their corresponding sequences of H 2n-1 following the above properties of the heights η(i).

This is done with the aid of a suitable tree of height 2n -1, in which each vertex is labeled by an integer of [n], so that: