Panagiotis Spentzouris

Iordanis Koutsopoulos

Grud Kasper

Tommy Vestergaard Madsen

Hansen

Kasper Grud

Tommy Vestergaard Hansen

Advertiser Bidding Prediction and Optimization in Online Advertising

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

Introduction

Online Advertising has evolved into a thriving business, with an average annual growth rate of 9.4%, and it is predicted to reach a total market size of $142.5 billion by 2021. In web-search advertising, advertisers cast bids in order to have their ads displayed in prominent positions in an ad list, next to organic results of a search-engine platform. Each ad contains a number of keywords. Advertisers set a budget they are willing to spend within a certain period of time, along with their bid for each keyword. Then, an auction is run by the platform and determines the ads to project and their ranking in the list, as well as relevant charges to advertisers, in terms of a cost-per-click (CPC). The form of auction that takes place is the Generalized Second-Price (GSP) one [START_REF] Varian | Online ad auctions[END_REF], and the rank of an ad is determined by the product of its bid and quality score, the latter being a cumulative estimate of the quality of the ad, including its design, text, graphics and landing pages.

Contrary to the vast amount of existing literature which considers the point of view of the advertising platform, we take the pragmatic viewpoint of an advertiser and study the problem of optimal bid selection across ads and across time, with the aim to maximize incoming traffic to the advertiser's page in terms of number of clicks, while we recognize that having data about other advertisers is practically infeasible. Maximization of inbound click traffic in the landing page is directly analogous to advertiser revenue, since a portion of the click traffic will end up completing a transaction and will produce revenue for the advertiser through a lead.

Each advertiser may leverage historical data about a number of quantities that are directly available to her, either through the advertising platform (e.g. Google AdWords), or through simple measurements at the advertiser side. These quantities are (i ), her bid values, i.e. how much the advertiser is willing to pay, (ii ) the associated CPC values, i.e. what the advertiser actually pays in the end of the auction, (iii ) the number of impressions of an ad, (iv ) the average position of an ad in the auction over a certain time interval, and (v ) the number of clicks received by the ad. The advertiser may then decipher the relationships between these measurable quantities and her bid. For example, the amount an advertiser actually pays for a certain ad position, i.e. the CPC value, depends on her bid, the quality score of the ad, as well as on the bids and quality scores of her competitors' ads. Also, the average position of an ad depends on its bid as well as the bids of competing ads. Finally, the number of clicks that an ad receives changes as a function of its average positions.

We capture the relations above by deriving regression models which may provide interesting insights, since they naturally contain in them hidden trends that would have been otherwise unavailable to the advertiser, such as the indirect impact of bidding behavior and quality scores of competing advertisers on the advertiser's ad position, CPC and number of clicks. The derived regression models provide the aforementioned relations in a simple analytic form that are then introduced to a convex optimization problem.

Our Contribution

Advertisers compete with each other for ad slots, over possibly several parallel auctions. Each advertiser has a set of ads, and each ad contains some keywords. For simplicity, we assume that each ad contains one keyword, and each keyword belongs to one ad. We study the problem of bid selection and allocation across ads and for different time intervals, in order to maximize the total incoming click traffic to the advertiser's landing page subject to a budget constraint. The contributions of our work to the literature are as follows.

-We consider the pragmatic viewpoint of a single advertiser that aims to best utilize data available to her in order to tune her bids for ads so as to maximize incoming click traffic to her landing page. -We follow a data-driven approach, where we build linear regression models about the relations of position versus bid, and CPC versus bid, and number of clicks versus position for each ad. These regression models demonstrate hidden trends about the competitor's behavior for different ads and different time intervals. -We incorporate the analytical expressions of the regression models into a convex optimization problem of bid selection across ads and across time, for maximizing the total number of ad clicks (and therefore the click traffic in the landing page), subject to a constraint on a maximum budget to be spent in a specific time period. The solution of the problem reveals interesting insights about the relation of the bid and the parameters of the regression models. -We validate our approach using real data. Our approach is shown to outperform other approaches that are either heuristic and data-agnostic or employ simple statistics on the data.

To the best of our knowledge, our work is among the first to follow the sequence of steps needed for an ad to be projected and possibly clicked through search advertising, and it uses machine-learning-based generated models to find hidden relations between the core ad attributes that are involved in these steps. We then incorporate them into a convex optimization problem the solution of which is the optimal bidding strategy for the advertiser of interest.

Model

Setup

We consider an advertising platform (e.g. Google AdWords) and a set of advertisers. We take the perspective of a single advertiser with a set A of N ads, where each ad has a quality score and contains some keywords. Without loss of generality, we assume that an ad includes exactly one keyword, so when we refer to an ad, we refer to its keyword. The advertiser has a budget B to be spent over Fig. 2. Our methodology in order to find the optimal bidding strategy in terms of maximizing total click traffic for the advertiser.

a period of time for bidding for her ads. The advertiser decides on bid vector b = (b 1 , . . . , b N ), where b i denotes the amount to bid for ad (keyword) i. The advertiser participates in different auction processes, one for each ad, and she competes with the same or similar keywords of other advertisers in order to have her ads displayed in as high a rank as possible in the lists. Each such auction takes as input the bids and quality scores of competing advertisers for the same or similar keywords and decides on the ranking of ads and the CPC to be paid by each ad. Finally, ads are displayed in the list, and they may be clicked by users who view the ad list next to organic search results. Once an ad is clicked, the user is taken to the landing page of the corresponding advertiser, from which a lead (i.e. a product purchase) might occur within some time interval. The advertiser thus earns a certain amount of revenue by each lead event.

For the advertiser of interest, let b i denote her bid for ad i. We implicitly assume that the auction for ad i is run once, hence b i is the amount that is bid for ad i. Let p i be the position of ad i in the ranked list. Although we know high positions in the ranked list have small absolute value, e.g. positions 1,2, we define the position in such a way that higher positions in the list are associated with high values. The position of each ad i in the list is equal to its negative, i.e. p i = -p i . Thus, if ad i is ranked first in the list, then p i = -1, if it is ranked second, then p i = -2, and so on. Also let c i be the cost-per-click (CPC) that is paid for ad i. Finally, let n i be the number of user clicks on ad i within some time interval.

There exist hidden dependencies between these quantities that we seek to capture in the sequel. First, the CPC value for ad i, c i , depends on bid value b i , i.e. the amount of money that the advertiser declares she is willing to pay for ad i. We denote this relation as c i = f i (b i ), where f i (•) denotes a continuous, and nondecreasing function. Along the same lines, position p i depends on the advertiser's bid b i , and let p i = g i (b i ) denote that relation, where g i (•) is a continuous, nondecreasing function. Finally, the number of clicks n i of ad i depends on its position, i.e. n i = h i (p i ) where h i (•) denotes again a non-decreasing function. 4

Dataset and linear regression models

The dataset that is readily available at the advertiser's side for each ad i is of the form:

D i = {(b d i , c d i , p d i , n d i ) : d = 1, . . . M i }, (1)
where M i is number of data points available in the dataset for ad i. For each ad i, we are interested in using dataset D i so as to build models for approximating the relations

f i (•), g i (•), h i (•).
We can use various machine-learning methods to derive models for the three relations above and then try to fit the models to the data in a way that reduces the total approximation error. Several methods are available, e.g. Neural networks, non-linear regression or linear regression models.

In this work, we adhere to linear regression models, because our objective is to demonstrate the advantages of fitting a data-driven model in an optimization problem and extract the benefits of optimization, rather than comparing different ML methods. Further, linear regression models provide a simple means to have the relations above in an analytic form so as to feed them in the optimization problem, to be presented in the next section, and derive interesting insights about the solution. Hence we consider the following three models:

c i = f i (b i ) = α i b i + β i , α i ≥ 0 (2)
p i = g i (b i ) = γ i b i + δ i , γ i ≥ 0 (3)
n i = h i (p i ) = λ i p i + µ i , λ i ≥ 0 (4)
where (α i , β i ), (γ i , δ i ) and (λ i , µ i ) are the parameters of the regression models to be computed from dataset D i for each ad i. The approach can be clearly extended in case the dataset is of different form.

The functions f i (•), g i (•), h i (•) above contain some hidden trends about the competing advertisers' bidding behavior as well. In the next section, we will include these models in an optimization problem that will give the optimal bid allocation policy in terms of total number of clicks for the ads of the advertiser.

Problem Statement and Formulation

Bidding for maximizing total number of clicks for ads

We are interested in finding the bid allocation policy b = (b 1 , b 2 , ..., b N ) that maximizes the advertiser's total click traffic for all ads. This is a key objective for an advertiser, since a percentage of this traffic will end up completing a transaction and produce revenue through leads. Presumably the advertiser participates in several auctions, one for each ad, and that auctions are independent from each other. We formulate the optimization problem as follows:

max N i=1 h i (p i ), (5)
subject to the constraint:

N i=1 h i (p i )f i (b i ) ≤ B, (6)
with b i ≥ 0 for i = 1, . . . , N , where B is a fixed amount of budget to be spent over a specific period of time.

Solution

Since both the objective function and the function involved in the constraint are increasing over the vector bid b, the constraint is satisfied with equality at the optimal solution b * . Thus, without loss of generality, it makes sense to consider the problem with an equality constraint. Let ω i = λ i γ i α i , and

ρ i = λ i (γ i β i + δ i α i ) + µ i α i .
This is a convex optimization problem, since both the objective function and the constraint are convex. Then we define the Lagrangian function, and after some algebraic manipulations we get the optimal solution as:

b * i = 1 2νa i - ρ i 2ω i (7)
Expression [START_REF] Zhang | Feedback Control of Real-Time Display Advertising[END_REF] shows that bid value b i decreases as a i (the slope of the straight line c i = f i (b i )) increases. Thus, if the CPC value of an ad increases with a high rate, then our algorithm will raise the bids of alternative ads whose CPC value grows at a slower rate. On the other hand, as the parameter ω i = λ i γ i α i increases, the bid value b i increases. Parameter ω i is the product of the slopes of the straight lines c i = f i (b i ), p i = g i (b i ), and n i = h i (p i ). If an ad earns higher positions and its number of clicks increases at a high rate, then our algorithm will "promote" this ad and increase its bid value. A second observation is that the product λ i γ i overpowers parameter α i , which is inversely analogous to the bid value b i . Therefore, the rate at which an ad's number of clicks increases or the rate at which it earns higher average positions weighs more in our algorithm than the rate at which the CPC value increases.

Remark 1: An alternative formulation could take into account the quality q i for each ad i ∈ A and include that as a factor in the objective function, which will now become N i=1 q i h i (p i ) and expresses the total weighted number of ads clicks, where the weight is the quality of each ad. By a similar reasoning as above, the optimal bid value for ad i is given as,

b * i = q i 2νa i - ρ i 2w i (8)
which tells us that the bid value of an ad grows proportionally to its quality score.

Remark 2: The optimization problem above considers budget allocation only across ads. We can have an enhanced formulation that would allow the advertiser to spend different amounts of budget at different time intervals. For example, it is natural to assume that the auction structure (i.e the number of competing advertisers and their bids) and its outcome will be different at different times of the day or at different days of the week (e.g. week-days or weekends) hence the bid would need to be adjusted as well. Let t = 1, . . . , T denote an index for different time intervals, where T is the total number for such time intervals. We will have another problem definition, where b i (t), p i (t), c i (t), n i (t) denote the bid, position, CPC and number of clicks of ad i at time interval t. Then, the budget allocation policy is given by vector (b(t) : t = 1, . . . , T ), with b(t) = (b 1 (t), . . . , b N (t)). Hence, the optimization problem becomes one of deciding on the bid policy that maximizes the total number of clicks across ads and across time, subject to a budget constraint.

Data Experiments

Dataset and Regression Models

We use a real dataset that was provided to us by an advertiser. The dataset consists of data from Google AdWords, and among other metrics it contains the bid values, CPC values, average position, and number of clicks for all keywords of the company. These metrics are available per day, for a time period of 6 months between August 2016 and February 2017. In our experiments we use 5 popular keywords/ads from the dataset since we are interested in demonstrating the benefits of the optimal policy. The properties of this policy are expected to hold for a greater number of keywords as well.

For the optimization problem ( 5)-( 6), we need 3 prediction models for each keyword and each time interval, i.e 30 regression models in total. The three regression models that we build for a keyword i and a time interval t = 1, 2 are, one for the estimation of the CPC c i = f i (b i ), one for the estimation of the average position p i = g i (b i ), and one for the estimation of the number of clicks

n i = h i (p i ).
In order to build the prediction models we used linear regression techniques from Machine Learning. Let us assume that our training dataset, for the regression model c = f (b) of a specific ad and a specific time interval is L = {b n , c n } M j=1 , where each b j (bid) is the data input and c j (CPC) the corresponding output, j is an index for each data point in the data set, and M is the size of the dataset, i.e. the number of different entries available to train each model. We will choose a polynomial of first degree as a model:

f (b j , w 0 , w 1 ) = w 0 + w 1 b j , (9)
and in order to find the parameters w 0 , w 1 that minimize the approximation error, we minimize over w the cost function of least squares. We measure the performance of the model by means of Root-mean-square-error. Finally, we use cross-validation in order to choose an appropriate value for the reguralization parameter.

Results

We run our experiment and compare the results for 4 different policies. We use the same set of 5 keywords for each policy, and initialize all bids at the value 0.1.

We run experiments for different values of budget B = {100, 200, 500, 1000, 1500, 2000}. The policies we study and compare are:

-Random bid allocation policy (Policy 1) is a baseline approach which allocates randomly a portion of the budget across ads, without taking into consideration any information about the performance of the ad. For this policy, we choose uniformly at random one of the five ads, and then we raise this ad's bid by 5%. We repeat this procedure until the portion of the budget is exhausted. We run experiments for portions equal to 10%, 20%, and 30% of the budget. For each portion we run the experiment 10 times, and keep as result the average number of clicks from these runs. -Inversely proportional to CPC allocation policy (Policy 2) allocates at each ad a portion of the budget that is inversely proportional to the average CPC value of ads, i.e. ads with high CPC get a smaller increase in their bids than ads with a lower CPC. Specifically, given the portion of the budget φ, which is the maximum amount of money that can be allocated across ads and across time without earning more clicks than the advertiser can pay, each ad's bid is raised by an amount 1/ci i 1/ci × φ. -Proportional to number of clicks allocation policy (Policy 3) allocates a portion of the budget φ, which is the maximum amount of money that can be allocated across ads and across time without earning more clicks than the advertiser can pay, that is proportional to the average number of clicks that ads get, i.e. ads which get a large amount of clicks have their bids increased by a bigger amount than ads which get less clicks. The increase of each ad's bid is ni i ni × φ. -Optimal policy (Policy 4) which is the outcome of the optimization problem ( 5), [START_REF] Xu | Lift-Based Bidding in Ad Selection[END_REF].

Comparison of Bid Allocation Policies. Figure 3 shows the total number of clicks to the advertiser's landing page under policy 1, for different budget amounts B = {100, 200, 500, 1000, 1500, 2000}, and different portions of the budget to be allocated across ads and across time equal to 10%, 20%, and 30%. The total number of clicks seems to be an increasing function over budget B for all cases except when the portion of budget is equal to 20%. Due to the randomness of the allocation, for this case and different budget amounts B = 1000, and B = 1500, it seems that the advertiser could earn more clicks with the smaller budget. Also, because of the random fashion in which we allocate a portion of the budget across ads and across time, it occurs that for a smaller portion of the budget we could have a better bid allocation, which could result in a larger amount of clicks to the advertiser's landing page.

In Figure 4, we depict the total number of clicks to the advertiser's landing page under three different policies. As expected, the total number of clicks seems to be an increasing function over the budget B for all policies. We see that by building regression models about the relations and by integrating them into the optimization problem ( 5), [START_REF] Xu | Lift-Based Bidding in Ad Selection[END_REF] we get the largest number of clicks, and thus larger revenue for the advertiser, compared to the other policies in the figure and the random policy in Figure 3. In Figure 5, we present the optimal bid selection for different amounts of budget B = {500, 2000}, different time intervals t = {1, 2}, i.e. week-day and week-end.

We observe that policy 2 tries to balance bids across all ads, and as a result it increases the bad performing ads' bids as well. We also see that under policy 3, there may be some cases of ads that are not currently getting many clicks but they have the potential to get more with less money than the current top performing ads. Nevertheless, this policy will ignore them because of the better performing in terms of clicks, ads. Despite this fact, policy 3 seems to perform better than policy 2. Also, we notice that policy 1 is performing similarly to the two other policies. Bid Allocation Policies weighted by Quality Score. In Fig. 6 we show the total number of clicks to the advertiser's landing page weighted by the quality score of each ad, under 4 different policies. The ads' quality scores are 9, 10, 8.8, 8, and 10 respectively. We see that the optimal policy outperforms the other 3 policies. Further, due to randomness, the random allocation can be benefited from increasing the bids of the ads with the biggest quality score values and thus outperform policy 3. Nevertheless, as the budget gets bigger policy 2 closes the gap in performance and finally outperforms policy 1 for a budget amount larger than 1700. This shows that for small amounts of budget, it is possible for a quality scoreagnostic policy to increase the bids of the non-top performing ads more than the bids of the top performing ones. If the budget is big enough though to increase the bids of all ads, all policies achieve better performance. Finally, policy 2 which tries to balance bids across all ads is again the worst performing policy.

Optimal bid allocation for maximum total ad quality. Figure 7, depicts the optimal bid selection for different budget values B = {500, 2000}, and different time intervals t = {1, 2}, i.e. week-day and week-end, for the objective that wants to maximize the total ad quality for the advertiser. We observe that for a small budget value B = 500 all the ads' bid values are increased compared to their values in Figure 5. However, for a bigger budget value B = 2000, the bid values of the ads with the best quality scores, i.e. the two ads with quality score q = 10, are significantly increased compared to their bid values in Figure 5, while the bid values of the rest are decreased. As expected, the optimal bidding policy for the objective that wants to maximize the total ad quality, decides to increase the bids of the ads that have the highest quality scores, while bid less for the rest. Furthermore, we observe that keyword 3 has a bid value b = 0 in both figures 5, 7. The main reason for that is the slope parameter λ of its regression model n = λp + µ, which show that keyword 3 gets clicks at a slower rate than the other keywords and thus our algorithm will not decide to increase its bid.

Related Work

Real-time bidding (RTB) represents the cutting-edge frontier of the computational advertising research, and an thriving research area in advertising together with the display ad network and search-based keyword advertising. In RTB, advertisers bid for an ad impression when a user visits a webpage. A repository with a nice taxonomy of recent works can be found at [START_REF]Paper Collection of Real-Time Bidding[END_REF].

Budget and bid allocation. In the context of display advertising, the authors in [START_REF] Cai | Real-Time Bidding by Reinforcement Learning in Display Advertising[END_REF], model the state transition via auction competition, and they build a Markov Decision Process framework for learning the optimal bidding policy to optimize the ads performance. In [START_REF] Zhang | Joint Optimization of Bid and Budget Allocation in Sponsored Search[END_REF], the authors define a revenue maximization problem, on an account level, by incorporating a probabilistic model to approximate the probability of winning a position given a price, and then they convert it into an integer optimization programming one. In work [START_REF] Li | Programmatic Buying Bidding Strategies with Win Rate and Winning Price Estimation in Real Time Mobile Advertising[END_REF] they study the problem of finding a bidding strategy in real time mobile advertising. First, they model the win rate using a logistic regression model, and then take the derivative of win rate estimation to generate the distribution of the winning price, and use the expected value of the distribution under the bid price as the winning price estimate. Then a bidding strategy is actually an optimization function that takes the input of expected revenue if winning the auction, win rate and winning price estimate, and generate the final bid price according to some pre-defined objective functions.

In another relevant work [START_REF] Zhang | Optimal Real-Time Bidding for Display Advertising[END_REF], the authors try to find the optimal bidding function that maximizes key performance indicator (KPI), i.e. the total number of clicks or revenue in Real Time Bidding (RTB) display advertising. They find a function that returns the probability of winning given a bid value, based on historic data, and then based on the form of the winning rate function they derive empirically a simple function that returns the bid value. Then, they feed these functions into an optimization problem that returns the optimal bid allocation for the advertiser. In a recent work [START_REF] Xu | Lift-Based Bidding in Ad Selection[END_REF], the authors try to use simple heuristic bidding policies to increase the number of clicks, and they set the bid for each ad impression proportionally to the increase of user's conversion rates.

Conclusion

We study the problem of bid selection for a single advertiser across ads, each of which is represented by linear regression models that are derived from real data. Our method showcases the benefit of feeding data-driven models into optimization problems, and to the best of our knowledge constitutes the major novelty of this paper.

Possible future steps in this work include studying the problem of maximizing advertiser's revenue across the whole path up to conversion, and pursuing sensitivity analysis to model and take into consideration the inaccuracy of the models.

Fig. 1 .

 1 Fig.1. A pictorial view of the sequence of steps in order for an ad to be projected and possibly clicked by users and produce revenue.

Fig. 3 .

 3 Fig. 3. Total number of clicks to the advertiser's landing page under random bid selection policy, for 3 different portions of the budget to be allocated across ads and time.

Fig. 4 .

 4 Fig. 4. Total number of clicks to the advertiser's landing page under 3 different bidding policies.

Fig. 5 .

 5 Fig. 5. Optimal bid selection for different amounts of budget B and different time intervals t.

Fig. 6 .

 6 Fig. 6. Total number of clicks weighted by the quality score of each ad under 4 different bidding policies.

Fig. 7 .

 7 Fig. 7. Optimal bid selection for different amounts of budget B and different time intervals t, for the objective that takes into account the quality score of ads.


Petia Koprinkova-Hristova

Miroslava Stefanova

Bilyana Genova
email: b.henova@abv.bg

Nadejda Bocheva

Echo State Network for Classification of Human Eye Movements during Decision Making

Keywords: Echo State Network, Time Series, Classification, Eye Tracking

The paper develops further a recently proposed author's approach for classification of dynamic data series using a class of Recurrent Neural Network (RNN) called Echo state network (ESN). It exploits the Intrinsic Plasticity (IP) tuning of ESN reservoir of neurons to fit their dynamics to the data fed into the reservoir input. A novel approach for ranking of a data base of dynamic data series into groups using the length of the multidimensional vector of reservoir state achieved after consecutive feeding of each time series into the ESN is proposed here. It is tested on eye tracker recordings of human eye movements during visual stimulation and decision making process. The preliminary results demonstrated the ability of the proposed technique to discriminate dynamic data series.

Introduction

In previous author's work [START_REF] Koprinkova-Hristova | Echo state networks for multi-dimensional dataclustering[END_REF] a novel approach for unsupervised clustering of static multidimensional data sets using a class of RNNs called Echo state networks [START_REF] Jaeger | Tutorial on training recurrent neural networks, covering BPPT, RTRL, EKF and the "echo state network" approach[END_REF][START_REF] Lukosevicius | Reservoir computing approaches to recurrent neural network training[END_REF] was proposed. Next it was successfully tested on numerous practical examples and the results were summarized in [START_REF] Koprinkova-Hristova | Multi-dimensional data clustering and visualization via Echo state networks[END_REF]. In [START_REF] Koprinkova-Hristova | Echo state networks in dynamic data clustering[END_REF] the approach was upgraded and applied for classification of dynamic data series too. The core of the approach, proposed first in [START_REF] Koprinkova-Hristova | Echo state networks for multi-dimensional dataclustering[END_REF], was to use the ESN to extract more informative features from multidimensional data sets. For this aim equilibrium states of the ESN reservoir neurons corresponding to every multidimensional data item presented to the ESN input were used. As it was shown in [START_REF] Koprinkova-Hristova | On effects of IP improvement of ESN reservoirs for reflecting of data structure[END_REF], the fitting of the ESN reservoir dynamics to reflect the input data structure can be achieved by an approach for ESN reservoir tuning called Intrinsic Plasticity (IP) [START_REF] Steil | Online reservoir adaptation by intrinsic plasticity for back-propagationdecoleration and echo state learning[END_REF][START_REF] Schrauwen | Improving reservoirs using intrinsic plasticity[END_REF] that is aimed at achieving the desired distribution of the ESN reservoir output.

Since the number of the new extracted features depends on the size of the ESN, question how to choose the most proper among them is still under investigation. Initially [START_REF] Koprinkova-Hristova | Echo state networks for multi-dimensional dataclustering[END_REF][START_REF] Koprinkova-Hristova | Multi-dimensional data clustering and visualization via Echo state networks[END_REF] it was proposed to choose only two of all possible neurons steady states based on their distribution in one-or two-dimensional space [START_REF] Koprinkova-Hristova | Multi-dimensional data clustering and visualization via Echo state networks[END_REF]. Next in [START_REF] Bozhkov | Reservoir computing for emotion valence discrimination from EEG signals[END_REF] we tried to extend the number of representative neurons until the accuracy of data clustering increases. However, this approach increased computational burden too much. That is why here another approach was proposed: to use the geometric size of the vector of all reservoir states and then to rank the time series data based on that single feature.

The approach was tested on dynamic series of eye movements data collected during psycho-physiological experiments with humans observing specific visual stimuli and making decisions. The preliminary results demonstrated the ability of the proposed approach to rank the time series of human eye movements in dependence on their characteristics.

The paper is organized as follows: next chapter describes briefly the ESN structure and its IP tuning and the newly proposed feature extraction approach; next the experimental set-up and the collected time series data are described; the results from classification by the proposed algorithm are presented and discussed in section four; the paper finished with concluding remarks and directions for future work.

2

Clustering Algorithm

Echo State Network and IP Tuning

ESN, shown on Fig. 1, is a type of recurrent neural network that belongs to the novel and fast developing family of reservoir computing approaches [START_REF] Jaeger | Tutorial on training recurrent neural networks, covering BPPT, RTRL, EKF and the "echo state network" approach[END_REF][START_REF] Lukosevicius | Reservoir computing approaches to recurrent neural network training[END_REF]. The ESN output for the current time instance k is the vector out(k) with size n out . It is a linear function f out (usually identity) of the vectors of the current states of the input in(k) (with size n in ) and the reservoir neurons X(k) (with size n X ):

              k X k in out W out f k y , (1)
Here W out is a trainable

) ( X n in n out n  
matrix. The neurons in the reservoir have a simple sigmoid output function f res (usually hyperbolic tangent) that depends on both the ESN input in(k) and the previous reservoir state X(k-1):

              1 k X res W k in in W res f k X (2)
Here W in and W res are The main reason for development of such a type of RNN is to simplify their training algorithm. However, it appeared that although non-trainable weights can be random, there is need to tune them initially. For this aim different approaches were proposed [START_REF] Jaeger | Tutorial on training recurrent neural networks, covering BPPT, RTRL, EKF and the "echo state network" approach[END_REF][START_REF] Lukosevicius | Reservoir computing approaches to recurrent neural network training[END_REF]. In [START_REF] Steil | Online reservoir adaptation by intrinsic plasticity for back-propagationdecoleration and echo state learning[END_REF][START_REF] Schrauwen | Improving reservoirs using intrinsic plasticity[END_REF] an algorithm called intrinsic plasticity (IP) was proposed. Its aim was to increase the entropy of the reservoir neurons outputs by minimization of the Kullback-Leibler divergence:

                      X d p X p r p X d p X p KL D log , (2)
that is a measure for the difference between the actual p(X) and the desired pd(X) probability distribution of reservoir neurons states X. It was proven that for the commonly used hyperbolic tangent at the reservoir neurons output the proper target distribution has to be the Gaussian one. For this aim two additional reservoir parameters, gain a and bias b (both vectors with size n X ), were introduced in [START_REF] Schrauwen | Improving reservoirs using intrinsic plasticity[END_REF] as follows:

                    b k X res W a diag k in in W a diag res f k X 1 (3)
The IP training is gradient descent algorithm [START_REF] Schrauwen | Improving reservoirs using intrinsic plasticity[END_REF] minimizing the Kullback-Leibler divergence by adjustment of the vectors a and b.

Classification Approach for Dynamic Data Series

In [START_REF] Koprinkova-Hristova | On effects of IP improvement of ESN reservoirs for reflecting of data structure[END_REF] was demonstrated that besides its initial aim, the IP tuning also fits the reservoir connections matrix to the structure of the input data presented to the ESN. Moreover, the equilibrium states of reservoir neurons corresponding to each one of the input data items used during IP tuning reflect the overall data structure [START_REF] Koprinkova-Hristova | Echo state networks for multi-dimensional dataclustering[END_REF][START_REF] Koprinkova-Hristova | Multi-dimensional data clustering and visualization via Echo state networks[END_REF]. Thus collected in this way features can be used for further classification or clustering.

In [START_REF] Koprinkova-Hristova | Echo state networks in dynamic data clustering[END_REF] it was demonstrated that the reservoir state X(N) reached after feeding of non-constant (time varying) sequence of inputs from in(0) to in(N) to the IP tuned ESN:

in(k) out(k)  reservoir W res X(k) W in W out                                                              b X res W a diag in in W a diag res f X b N X res W a diag N in in W a diag res f N X b N X res W a diag N in in W a diag res f N X 0 0 1 2 1 1 1  (4)
depends on dynamic characteristics of the time series in and can be exploited as a set of classification features.

Since the choice of neurons whose states are the best feature set for each particular data subject to classification or clustering is non-trivial [START_REF] Koprinkova-Hristova | Multi-dimensional data clustering and visualization via Echo state networks[END_REF][START_REF] Bozhkov | Reservoir computing for emotion valence discrimination from EEG signals[END_REF], here we propose another approach: to calculate the size of the vector containing all collected reservoir neurons states:

            N x N x N x N X X n i N i x R X n  2 1 , 1 2     (5)
and to use it as single discriminating data feature.

Experimental Set-up

The time series data used to test the idea described above were collected by eye tracing device that recorded the human eye movements during a behavioral experiment performed with the participation of volunteer human subjects observing series of visual stimuli.

Each stimulus is composed by a sequence of consecutive frames. A frame contains of 50 dots presented in a circular aperture with a radius of 7.5cm in the middle of the computer screen. The dots were grouped in 25 pairs placed at 2cm distance from each other. Each pair of dots had a limited lifetime of 3 frames. On every frame one-third of the pairs changed position. Each frame lasted 33ms. The orientation of the virtual lines connecting the dots in 18 pairs intersected in a common point considered as the center of each frame, while the rest 7 pairs had random orientations. The mean position of the centers of all frames in a stimulus sequence determines its "imaginary" center. We generated 14 different types of stimuli having centers at 7 positions shifted left and 7 positions shifted right from the screen midpoint. All shifts were in horizontal direction and varied between 0.67cm and 4.67cm with step of 0.67cm. Ten different patterns for each center position were generated.

The stimuli were presented on a gray screen with mean luminance 50cd/m 2 using 20.1′′ NEC MultiSync LCD monitor with NvidiaQuadro 900XGL graphic board at a refresh rate of 60 Hz and screen resolution 1280×1024 pixels. The experiments were controlled by a custom program developed under Visual C++ and OpenGl.

The subject sat at 57cm from the monitor screen. Each stimulus presentation was preceded by a warning sound signal. A red fixation point with size of 0.8cm appeared in the center of the screen for 500ms. The stimuli were presented immediately after the disappearance of the fixation point. The Subject's task was to continue looking at the position where the fixation point was presented until he/she made a decision where the center of the pattern was and to indicate this position by a saccade (fast eye movement). The subjects also had to press the left or the right mouse button depending on the perceived position of the center -to the left or to the right from the middle of the screen. If the subject could not make a decision during the stimulus presentation (3.3s for 100 consecutive frames), the stimulus disappeared and the screen remained gray until the subject made a response.

The eye movements of the participants in the experiment were recorded by a specialized hardware -Jazz novo eye tracking system (Ober Consuting Sp. Z o.o.). All recordings from all the sensors of the device for one session per person were collected with 1 KHz frequency and the information is stored in files. These include: the calibration information; records of horizontal and vertical eye positions in degrees of visual angle eyex and eyey; screen sensor signal for presence/absence of a stimulus on the monitor; microphone signal recording sounds during the experiment; information about tested subjects (code) and type of the experimental trail for each particular record.

The raw data were processed to extract only the records during presence of a stimulus on the screen. The data between the stimuli was excluded since it is not relevant to the eye movements during task performance.

Three age groups took part in the experiment: young (from 20 to 35 years), elderly (from 57 to 84 years) and middle age group (from 25 to 55 years). From all collected experimental data we observed big variety of eye movement behaviors varying not only between three age groups but also within each group. So it was very hard to classify test subjects only on the basis of this information. Hence we decided to try whether the proposed above dynamic data discrimination approach can yield some reasonable results.

The input to the ESN was two dimensional vector composed by the visual angles data series recorded during presence of a stimulus on the screen, i.e.

        k eye k eye k in y x



. We tuned three ESNs with reservoir sizes 10, 50 and 100 neurons using the IP algorithm described above. The extracted in this way feature of each dynamic data series was R calculated according to the equation ( 5). Thus the tested subjects were ordered based on the obtained value of R from their recorded eye movement data series.

Classification Results and Discussion

First we selected a representative group of four experienced test subjects from different age groups. These subjects took part in experimental set-up preparation so they were able to perform the behavioral tasks strictly and their eye movement recordings were clear from outliers due to improper behavior like looking to the mouse before clicking when decision was taken or keeping fixation. Such noisy behavior was observed with other volunteer subjects especially during their first trails.

Fig. 2 represents the eye movements' data series collected from these four "experienced" test subjects who performed first the described above experiment. Subject 1 is the youngest, subject 3 is middle aged and other two subjects (2 and 4) belong to the elderly group. We can easily distinguish the middle aged subject while similarities between subjects 2 and 4 are not so obvious. From Fig. 3, representing the variances of the data series, we can conclude that age differentiation by this characteristic is also a hard task even for such small group. Then we applied described above classification approach to these data series. In order to prove the expected effect of IP tuning of the ESN reservoir, we compare the obtained feature value before (Fig. 4) and after (Fig. 5) its application to the three randomly generated ESN reservoirs containing 10, 50 and 100 neurons.

From Fig. 5 we can conclude that IP tuning definitely helps to classify our four subjects by their age no matter of the reservoir size. group. This was confirmed by our algorithm according to which the middle aged subject 3 is clearly differentiated from the other three experienced subjects of different ages. Moreover, the two elderly subjects 2 and 4 were classified as close to each other while the younger person 1 was clearly differentiated too. Then we proceed with data recorded from 18 volunteers from three age groups as shown on figures 6, 7 and 8 respectively. Fig. 9 represents the variances of all data series from figures 6-8.

Although we observe some similarities and differences between the three groups, there are also significant dissimilarities between subjects from the same group. For example, we can conclude that young subject 6 behaved significantly different from other members of Group 1; in Group 2 we can distinguish subjects 10 and 11 from others while in Group 3 subject 15 seems to have different eye movement behavior.

Here we decided to IP tune two ESN reservoirs (denoted here as ESN1 and ESN2) from each size since initial connections of ESN were randomly generated, and to compare the obtained results. Fig. 10 represents the achieved values of R from ESN1 and ESN2 as well as its mean value. It is clear that for both initial ESN reservoirs we obtained similar results after IP tuning, especially for the middle group.

Fig. 11 represents the corresponding rank number of the subjects in dependence on the value of R from Fig. 10. The obtained subject order is similar for both ESN initial reservoirs. For most of the subjects two generations of ESN yielded the same rank. The same are the results from the ranking according to the mean value of R. The bigger differences are observed for the first and second subjects only.

However the obtained results showed that differentiation of tested subjects by age using only recorded eye movement behaviors is not possible since there are subjects from different age groups that received close ranks. Nevertheless, the approach seems promising for classification of types of eye movements during decision making that could be related to other psycho-physiological peculiarities of the tested subjects.

Conclusions

In conclusion, the proposed approach for classification of time series using the geometric size of the ESN reservoir vector state appears promising since it decreases significantly computational burden of the algorithm. Concerning the classification of human subjects based on their eye movements, it is obvious that obtained ranking is not enough to separate the people into the groups related to their age. Since the achieved ranking of subjects remain stable using different initializations and sizes of ESN reservoir, we can search other similarities between tested subjects that this classification might reveal.

Another explanation of the reported results might be presence of outliers in the collected data since in these preliminary investigations the raw data was used. Hence further refinement of the collected experimental data base could help to reveal some age-related similarities of the recorded eye movements.

Another direction for future work can be inclusion of the additional information collected during the experiments like amplitude, velocity and acceleration of saccades performed during decision making as well as accuracy of the persons' response and the corresponding reaction time for each individual stimulus. All these characteristics can serve as features to support classification of test subjects in our future work.

Fig. 1 .

 1 Fig. 1. Echo state network structure.

Fig. 2 .

 2 Fig. 2. Experienced 4 subjects' recordings from eye tracker.

Fig. 3 .

 3 Fig. 3. Variances of eye movement coordinates from Fig. 2.

Fig. 4 .

 4 Fig. 4. Size of the reservoir state vector R achieved after presentation of input time series for four experienced test subjects before IP tuning.

Fig. 5 .

 5 Fig. 5. Size of the reservoir state vector R achieved after presentation of input time series for four experienced test subjects after IP tuning.

Fig. 6 .

 6 Fig. 6. Eye movements recorded from young test subjects.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Eye movements recorded from the middle aged test subjects.

Fig. 9 .

 9 Fig. 9. Variances of eye movement coordinates for all test subjects.

Fig. 10 .

 10 Fig. 10. Size of the reservoir state vector R for the three ESN reservoirs and all test subjects. ESN1 and ESN2 denote the results from the first and second ESN respectively; mean is the mean value between ESN1 and ESN2.

Fig. 11 .

 11 Fig. 11. Rank of the subjects according to the size of the reservoir vector R for the three ESN reservoirs and all test subjects. ESN1 and ESN2 denote the results from the first and second ESN respectively; mean is their mean value.

R ESN reservoir size = 100 ESN 1 ESN 2 mean ESN 1 ESN 2 mean ESN 1 ESN 2 mean

											ESN reservoir size = 10
			3
	R	2.5
			2
			1 1.5	2	3	4	5		6	7	8	9	10	11	12	13	14	15	16	17	18
										ESN reservoir size = 50
			5
	R	4
			3
			1 2	2	3	4	5	6	7	8	9		10	11	12	13	14	15	16	17	18
			8
			6
			4
			1 2	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
											test subject number
											ESN reservoir size=10
	No	0 2 4 6 8 10 18 16 14 12	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
											ESN reservoir size=50
	No	0 2 4 6 8 10 18 16 14 12	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18
	0 2 4 6 8 10 12 14 16 18	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18

test subject number No ESN reservoir size=100 ESN 1 ESN 2 mean ESN 1 ESN 2 mean ESN 1 ESN 2 mean

Acknowledgment

The reported work is a part of and was supported by the project № DN02/3/2016 "Modelling of voluntary saccadic eye movements during decision making" funded by the Bulgarian Science Fund.


Amelia Bȃdicȃ
email: ameliabd@yahoo.com

Costin Bȃdicȃ

Daniela Dȃnciulescu

Doina Logofȃtu
email: logofatu@fb2.fra-uas.de

Greedy Heuristics for Automatic Synthesis of Efficient Block-Structured Scheduling Processes from Declarative Specifications

Keywords: Greedy algorithm, process model, ordering constraints, optimization

This paper introduces a new Greedy heuristic algorithm for the automatic synthesis of block-structured scheduling processes that satisfy a given set of declarative ordering constraints, as well as basic theoretical results that support the correctness of this algorithm. We propose two heuristics that can be used with this algorithm: hierarchical decomposition heuristic and critical path heuristic. We also present initial experimental results supporting the effectiveness and efficiency of our proposed algorithm and heuristics.

Introduction

There are many formalisms for the specification of business process models. Blockstructured models have certain advantages compared with other approaches [START_REF] Leemans | Discovering Block-Structured Process Models from Event Logs -A Constructive Approach[END_REF].

It is useful and quite intuitive to declaratively specify desired properties of process models. We are interested in constructing process models that are consistent with the given declarative specification [6]. This problem has practical applications in scheduling tasks encountered in manufacturing systems [START_REF] Mrasek | Automatic Generation of Optimized Process Models from Declarative Specifications[END_REF].

Manual construction of large process models satisfying a set of ordering constraints is almost impossible or at least not scalable. Automatic generation based on exhaustive exploration of the space of possibilities is difficult because of the huge number of potential candidates. The only feasible solution is to design automatic approaches based on efficient heuristic algorithms that are able to drastically prune the huge search space.

In this paper we focus on scheduling processes. In this case, the declarative specification defines the scheduling constraints. We are interested in determining optimal or at least, as efficient as possible, block-structured scheduling processes that satisfy the scheduling constraints. The optimization criterion requires the minimization of the total completion time. Optionally, we can add other constraints, like for example imposing upper bounds for the amount of parallel work. This constraint may result from the practical restriction regarding the limited availability of certain resources. In particular: i) our processes are defined only using sequential and parallel composition; ii) each activity must have exactly one instance in the schedule.

Our work was mainly influenced by previous results of [START_REF] Mrasek | Automatic Generation of Optimized Process Models from Declarative Specifications[END_REF]5]. Nevertheless, our results are different in many aspects. Most important, our heuristics are deterministic and different. We are using the hierarchical decomposition of a graph, while [START_REF] Mrasek | Automatic Generation of Optimized Process Models from Declarative Specifications[END_REF]5] are based on the more complex modular decomposition. We also provide theoretical results to support our work. Finally, we performed experiments with larger graphs, and our preliminary results suggest that our algorithm might be faster.

Note that there are many theoretical studies on evolutionary algorithms and randomized (meta/hyper) heuristics applied to combinatorial optimization algorithms [START_REF] Demertzis | Adaptive Elitist Differential Evolution Extreme Learning Machines on Big Data: Intelligent Recognition of Invasive Species[END_REF]. Such works could be considered for the further expansion of our results by comparison of our method with different, but related approaches.

Process Models

Let us consider a finite nonempty set of activities Σ. A trace t ∈ Σ * is a sequence of zero or more activities 3 . The length of a trace t = a 1 a 2 . . . a n is n and this is denoted as |t| = n. The empty trace is denoted by ε and |ε| = 0. For each nonempty trace t = a 1 a 2 . . . a n we define: i) the head of t as head(t) = a 1 , and ii) the tail of t as tail(t) = a 2 . . . a n .

A language L ⊆ 2 Σ * is defined as a set of traces. We can define certain operations with languages. The sequential composition of two languages L 1 and L 2 denoted by L 1 → L 2 , is defined as follows:

L 1 → L 2 = {w = l 1 l 2 | l 1 ∈ L 1 and l 2 ∈ L 2 }
This notation can be extended for a trace t and a language L as: t → L = {t} → L. The parallel composition of two traces t 1 and t 2 , denoted by t 1 t 2 , is defined as:

-For each nonempty trace t we have: t ε = ε t = {t} -For each nonempty traces t 1 and t 2 we have:

t 1 → t 2 = (head(t 1 ) → (tail(t 1 ) t 2 )) ∪ (head(t 2 ) → (t 1 tail(t 2 )))
The parallel composition L 1 L 2 of two languages L 1 and L 2 is now defined as:

L 1 L 2 = ∪ t 1 ∈L 1 ,t 2 ∈L 2 t 1 t 2
Let us consider the set {→, |, } of three binary operators used for constructing blockstructured processes. The operator → denotes sequential composition, the operator | denotes the nondeterministic choice, and the operator denotes parallel composition.

Let us denote with a, b, c, . . . the activities of Σ and with P, Q, R, . . . process terms. Process terms can be defined recursively as follows:

P ::= a | P → Q | P | Q | P Q
The language L(P) of process P is recursively defined as follows:

-L(a) = {a} -L(P → Q) = L(P) → L(Q) -L(P | Q) = L(P) ∪ L(Q) -L(P Q) = L(P) L(Q)
Operator has higher precedence, operator → has middle precedence, and operator | has lower precedence. All operators are associative, while and | are also commutative.

Process terms represent models of processes and they can be graphically depicted as trees or as block-structured flowcharts, as shown in Figure 1. In what follows we focus on process models with the following particularities:

-They represent sets of possible activity schedules. A schedule must contain exactly one instance of each activity. -They use sequential (→) and parallel ( ) operators. Scheduling processes are deterministic, explaining why nondeterministic choice is not used in their definition.

Rigorously defining scheduling processes requires the introduction of the support set supp(P) of a process P that denotes the set of activities that occur in process P.

A block-structured scheduling process is recursively defined as follows:

-If a is an activity then a is also a process such that supp(a) = {a}.

-If P and Q are processes such that supp(P) ∩ supp(Q) = ∅ then P → Q and P Q are processes with supp(P → Q) = supp(P Q) = supp(P) ∪ supp(Q).

For example, processes a c → b and a (c → b) are well-formed, and:

-supp(a c → b) = supp(a (c → b)) = {a, b, c} -L(a c → b) = {acb, cab} -L(a (c → b)) = {acb, cab, cba}
It is not difficult to observe that if P is a well-formed block-structured scheduling process then all its traces t ∈ L(P) have the same length |t| = |supp(P)|.

3 Declarative Specification of Ordering Constraints

Activity Ordering Graph

Based on domain-specific semantics, one can impose ordering constraints of the activities of a process. For example if two activities are independent and there are enough resources to be allocated to each of them then those activities can be scheduled for parallel execution. However, if an activity depends on the output produced by another activity, then the first activity can be scheduled for execution only after the completion of the second activity, i.e. there is a sequencing constraint between their execution order. Finally, if two activities define distinct action options then their execution is incompatible, so it cannot occur within the same schedule, i.e. they are mutually exclusive.

The ordering constraints imposed on each trace of a scheduling process are declaratively specified using an activity ordering graph G = V, E [5] such that:

-V is the set of nodes and each node represents an activity.

-E ⊆ V × V is the set of edges. Each edge represents an ordering constraint. Set E is partitioned into two disjoint sets E → and E with the following meaning:

• Set E → specifies sequential ordering constraints. If (u, v) ∈ E → then activity v cannot occur in a schedule without being preceded by activity u. E → is a partial ordering, i.e. it is transitive and antisymmetric, so it cannot define cycles.

• Set E specifies mutual exclusion constraints. If (u, v) ∈ E then activities u and v are incompatible, so they cannot occur within the same schedule. Set E defines a symmetric relation.

Intuitively, satisfaction of mutual exclusion constraints requires the availability of nondeterministic choice operator in process definition. As we assumed that this operator is not available for scheduling processes, we will now focus only on sequential ordering constraints, i.e. we assume that E = ∅ so E = E → . This means that the ordering graph is a directed acyclic graph with arcs defining sequential ordering constraints.

If t = a 1 a 2 . . . a n is a trace of a scheduling process and u, v are two activities of t then u precedes v in t, i.e. u t → v if there are 1 ≤ i < j ≤ n such that a i = u and a j = v. Let G = V, E be an ordering graph and let t be a trace containing all the activities of V with no repetition. Then t satisfies G, written as t | = G, if and only if E → ⊆ t →. This means that trace t cannot contain activities ordered differently than as specified by G.

The language L(G) of an ordering graph G is the set of all traces that satisfy G, i.e:

L(G) = {t | t | = G}
Let P be a scheduling process and let G = V, E be an ordering graph. P satisfies G written as P | = G, if and only if:

-L(P) ⊆ L(G), i.e. each trace of P satisfies G, and supp(P) = V, i.e. all the activities of V are relevant and occur in P.

The set of processes P such that P | = G is nonempty, as it contains at least one sequential process defined by the topological sorting of G.

Optimal Scheduling Processes

Each activity has an estimated duration of execution that is represented using a function d : Σ → R + . The duration of execution d(P) of a process P is defined as follows:

-

If P = a then d(P) = d(a). -d(P → Q) = d(P) + d(Q). -d(P Q) = max {d(P), d(Q)}.
The minimum duration of execution of a process that satisfies a given ordering graph G, denoted with d MIN (G), is defined as:

d MIN (G) = min P| =G {d(P)}
An optimal scheduling process that satisfies a given ordering graph G is a process P * with a minimum duration of execution, i.e. it satisfies:

-P * | = G, and -d(P * ) = d MIN (G).
There is a finite and nonempty set of processes that satisfy an ordering graph G, so the optimal scheduling process trivially exists. Moreover, as there is an exponential number of candidate processes satisfying G, we postulate that the computation of the optimal scheduling process is generally an intractable problem. Therefore, we will be focusing on developing efficient heuristic algorithms that are able to produce "suboptimal" or "good enough" scheduling processes using a reasonable computational effort.

Heuristics for Suboptimal Processes

We introduce two heuristics that are used to derive an efficient Greedy heuristic algorithm for computing a suboptimal scheduling process satisfying an ordering graph.

Hierarchical Decomposition Heuristic

Let G = V, E be an ordering graph. Remember that G is a directed acyclic graph defining the sequential ordering constraints imposed on a scheduling process.

-For each node v ∈ V we define the set I(v) of input neighbors of v as follows:

I(v) = {u ∈ V | (u, v) ∈ E}.
-For each node v ∈ V we define the level l(v) of v as a function l : v → N such that:

• If I(v) = ∅ then l(v) = 0. • If I(v) ∅ then l(v) = 1 + max u∈I(v) {l(u)}. -The height l(G) of graph G is defined as l(G) = max v∈V {l(v)}.
-

If m = l(G) ≥ 0 then the family of m + 1 sets {V 0 , V 1 , . . . , V m } defined as V i = {v | l(v) = i} for all 0 ≤ i ≤ m is a partition of V. If G i is the subgraph of G in- duced by V i then the family of graphs {G 0 , G 1 , . . . , G m } is known as the hierarchical decomposition of G.
Proposition 1. (Hierarchical Decomposition Process) Let G = V, E be an ordering graph. The hierarchical decomposition process P HD (G) associated to G is defined as:

-P i = v∈V i v for all 0 ≤ i ≤ m. -P HD (G) = P 0 → P 1 → • • • → P m .
Then Figure 2 shows an ordering graph G 1 , and two processes P 1 and P 2 such that G 1 | = P 1 and G 1 | = P 2 . The hierarchical decomposition of G 1 is induced by the partition of its vertices {{a, c}, {b}}, so we can easily notice that P 1 is the hierarchical decomposition process of G 1 . Observe that:

-d HD (G 1 ) = d(P 1 ) = max {d(a), d(c)} + d(b) -d(P 2 ) =

Critical Path Heuristic

Observe that an activity u cannot start unless all the neighboring activities from the input set I(u) are finished. This time point is denoted with start(u). Activity u that started at start(u) will finish at time f inish(u) = start(u)+d(u). The values start(u) and f inish(u) for each activity u ∈ V can be computed using the critical path method [2], as follows:

-If I(u) = ∅ then start(u) = 0 and f inish(u) = d(u). -If I(u) ∅ then start(u) = max v∈I(u) { f inish(v)} and f inish(u) = start(u) + d(u).
The maximum value of the finishing time of each activity, known as critical path length, is a lower bound for the duration of execution of the optimal scheduling process.

Reducing the Duration of Execution

Analyzing Figure 2, we can observe that the duration of execution of the hierarchical decomposition process can be reduced by doing a transformation that pushes the parallel composition operations upper in the process tree. However, this transformation is not always possible. We now provide sufficient conditions that enable the transformation and guarantee that the duration of execution of the resulted process is lower than of the original process. Referring at Figure 2, the key observation is that the set of nodes of graph G 1 can be partitioned in two subsets U 0 = {a, b} and U 1 = {c} such that there are no arcs cross-linking nodes in U 0 to nodes in U 1 or nodes in U 1 to nodes in U 0 . Note that such a decomposition is not possible for the graph G 2 from Figure 2. We consider the most general situation of reducing the duration of execution of a process (P 1 P 2 ) → (Q 1 Q 2 ). Similar results can be obtained for the processes of the form (

P 1 P 2 ) → Q and P → (Q 1 Q 2 ). Proposition 3. (Reducing the Duration of Execution) Let P = (P 1 P 2 ) → (Q 1 Q 2 )
and let G = V, E be an ordering graph such that P | = G. Let us also assume that

((supp(P 1 → Q 1 ) supp(P 2 → Q 2 )) ∩ E = ∅ 4 .
Then it follows that:

-Process P = (P 1 → Q 1 ) (P 2 → Q 2 ) is well-formed, -P | = G,

Automatic Synthesis Algorithm

Let G = V, E be an ordering graph and let U be the undirected graph obtained by removing the orientation of arcs of graph G. We denote with G(W) and U(W) the subgraphs of G and U induced by a subset W ⊆ V of nodes.

Let {V 0 , V 1 , . . . , V m } be the partition of node set V defined by the hierarchical decomposition of G. We define the following sets of nodes:

-W 0 = V 0 -W 1 = W 0 ∪ V 1 . . . -W m = W m-1 ∪ V m = V Let U i = U(W i ) and G i = G(W i ) for each 0 ≤ i ≤ m.
Each undirected graph U i can be partitioned into connected components that induce the partition {U 1 , U 2 , . . . , U k i } of the set W i of nodes such that k i > 1. This situation is intuitively described in Figure 5.

Following the result of Proposition 3, the hierarchical decomposition process P defined for subgraph G i can be transformed into process P such that:

-P = k i j=1 P j -supp(P j ) = U j for all 1 ≤ j ≤ k i -d(P) ≥ d(P )
Consider for example the sample ordering graph G 3 from Figure 6. The partition of nodes corresponding to the hierarchical decomposition of

G 3 is {V 0 , V 1 , V 2 } = {{a, c}, {b, d}, {e}} and its height is m = 2. The hierarchical decomposition process of G 3 is P 6 = (a c) → (b d) → e.
Its duration of execution is 46.

We observe that for i = 1 the set

W 1 = V 0 ∪ V 1 = {a, b, c
, d} can be partitioned into {{a, b}, {c, d}}, so k 1 = 2. Using this observation we determine the transformed process It follows that by applying our proposed transformation we were able to significantly reduce the duration of execution of process P 6 from d(P 6 ) = 46 to d(P 7 ) = 36. We can combine this transformation with the hierarchical decomposition heuristic d HD provided by Proposition 1 or with the critical path heuristic d CP provided by Proposition 2 to design an efficient Greedy algorithm for the automatic synthesis of a suboptimal scheduling process that is consistent with a declarative specification.

V 0 V 1 . . . V i W i ⇒ U k i . . . U 2 U 1 W i
Let G = V, E be an ordering graph. The algorithm can be defined as a function proc(W, G(W)) that takes a subset of nodes W ⊆ V, the subgraph G(W) of G induced by W and returns a suboptimal process that satisfies G(W).

Let {V 0 , V 1 , . . . , V m } be the partition of node set V defined by the hierarchical decomposition of G. Function proc(V, G(V)) is recursively defined as follows:

-If m = 0 then proc(V, G(V)) = v∈V v. -If m > 0 and V 0 = {v} is a singleton set then proc(V, G(V)) = v → proc(V \ {v}, G(V \ {v})).
-If m > 0 and V 0 has at least two elements then for each 0 ≤ i ≤ m determine the number k i of the sets of the partition of set W i induced by the connected components of the undirected graph U i obtained from the directed graph G i . We have

k 0 ≥ k 1 ≥ • • • ≥ k m ≥ 1.
Let i be the largest index for which k i > 1. Such an index always exists as k 0 = |V 0 | > 1. Select an index 0 ≤ j ≤ i for which the estimated duration of execution of the "synthesized process" (to be defined in what follows) is minimized.

We now recursively define the "synthesized process" and its estimated duration of execution, in terms of function proc. Let G = V, E be an ordering graph, let {V 0 , V 1 , . . . , V m } be the partition of node set V defined by the hierarchical decomposition of G, and let us assume that m > 0 and |V 0 | > 1. The "synthesized process" P j and its estimated duration of execution d G-ES T (P j ) with ES T ∈ {HD, CP} is:

-If j = 0 then P 0 = v∈V 0 → proc(V \ V 0 , G(V \ V 0 )) and is duration of execution is estimated to d G-ES T (P 0 ) = max{d(v)} v∈V 0 + d G-ES T (G(V \ V 0 )).
-If 0 < j < m then let us consider the partition {Y 1 , Y 2 , . . . , Y k j } of W j . Then P j = (

k j i=1 proc(Y i , G(Y i ))) → proc(V \ W j , G(V \ W j )) and is duration of execution is estimated to d G-ES T (P j ) = max k j i=1 {d G-ES T (Y i )} + d G-ES T (V \ W j ). -If j = m then let us consider the partition {Y 1 , Y 2 , . . . , Y k m } of W m = V. Then P m = k m i=1 proc(Y i , G(Y i ))
(G) ≥ d G-ES T (G) ≥ d MIN (G) ≥ d CP (G).

Experimental Evaluation

We implemented our algorithm in Standard C using the 64-bit GCC compiler, version 5.1.0 and tested it on a x64-based PC with Intel(R) Core(TM) i7-5500U CPU at 2.40GHz running Windows 10. In this section we present the experimental results that we obtained with this implementation. The experiment was organized as follows:

-We randomly generated a number of directed acyclic graphs of increasing sizes representing ordering constraints, as well as random durations of execution for each activity of the graph. The parameters of a data set are: number n of graph nodes, number ng of generated graphs, minimum and maximum durations dmin and dmax of each activity, and the density factor f ∈ [0, 1] of the graph. The higher is this factor the more dense is the graph. Value of f is given as a percentage. -For each graph G we estimated the basic metrics given by the hierarchical decomposition heuristic d HD (G) and by the critical path heuristic d CP (G). -For each graph G we computed the suboptimal scheduling process that satisfies G using the Greedy heuristic algorithm proposed in Section 4.4, in two variants: using the hierarchical decomposition heuristic and respectively using the critical path heuristic, to confirm the result claimed by Proposition 4, and to compare the results obtained for d G-HD and d G-CP .

The graph data sets were generated for the following values of the parameters: ng = 100, n ∈ {10, 50, 150, 300, 500, 700}, dmin = 1, dmax = 20, and density factor f ∈ {15%, 30%, 45%, 60%, 75%}. For each test we recorded the total execution time and the values of the metrics of interest. We labelled each data set to reflect its number of nodes and density. For example if n = 500 and f = 30% then the label is 500-30. 1 presents the total execution time of running the synthesis algorithm for each data set. We observe that increasing the number of nodes, as well as the density, determines the increase of the execution time. Note that these times cover the processing of batches of 100 graphs. This means for example that the average time to process one graph of the 700-75 data set is approximately 1 second, i.e. our algorithm is quite fast. Secondly, that results of both experiments show that CP heuristic performs better than HD heuristic for almost all the graphs of the data set (there are few exceptions difficult to observe on the figures). Thirdly, the heuristics CP and HD tend to give closer results for higher density ordering graphs, as can be noticed by comparing the "closeness" of the cost values obtained for G-CP and G-HD for each data set 700-30 and 700-60.

Conclusions

We proposed a new Greedy algorithm for the automatic synthesis of block structured scheduling processes that satisfy given declarative ordering constraints. We presented basic theoretical results that support the correctness of this algorithm. We proposed two heuristics that can be used with this algorithm: hierarchical decomposition and critical path. Our initial experimental results support the effectiveness of our proposals and suggest that the critical path heuristic performs better.

Fig. 1 .

 1 Fig. 1. Tree model of process a → b c (left) and its equivalent block-structured model (right)

Fig. 2 .

 2 Fig.2. Ordering graph G 1 (left), process P 1 (middle) and process P 2 (right)

  max {d(a) + d(b), d(c)} Clearly d(P 1 ) ≥ d(P 2 ) and P 2 is optimal (other satisfying processes are strictly sequential, incurring a higher duration of execution). But note that if d(a) ≥ d(c) then d(P 1 ) = d(P 2 ) = d(a) + d(b) so the optimal scheduling process has duration d HD (G 1 ) which shows that we can have equality in the inequality resulted from Proposition 1. However, if d(a) < d(c) the optimal scheduling process has duration d(P 2 ) = max {d(a) + d(b), d(c)} < d HD (G 1 ) = d(c) + d(b).

Proposition 2 .

 2 (Critical Path) Let G = V, E be an ordering graph and let d CP (G) be its critical path length. Then d CP (G) is a lower bound of the duration of execution of the optimal scheduling process d MIN (G), i.e. d MIN (G) ≥ d CP (G).

Fig. 3 .Figure 3 Fig. 4 .

 334 Fig. 3. Ordering graph G 2 (left) and process P 3 (right)

  and d(P) ≥ d(P ).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Transformation to reduce duration of execution.

  and is duration of execution is estimated tod G-ES T (P m ) = max k m i=1 {d G-ES T (Y i )}. Proposition 4. (Duration of Execution of Greedy Suboptimal Processes) Let d G-ES T (G)be the duration of execution of the suboptimal process that was computed with the Greedy algorithm using heuristic ES T ∈ {HD, CP}. Then this process satisfies ordering graph G and d HD

Fig. 7 .

 7 Fig. 7. Comparison of costs for the 700-30 data set

Fig. 8 .

 8 Fig. 8. Comparison of costs for the 700-60 data set

Table 1 .

 1 Total execution time in seconds for processing each data set Table

	# nodes / density 15% 30% 45% 60%	75%
	10	0.019 0.009 0.012 0.010 0.012
	50	0.191 0.202 0.186 0.208 0.198
	150	1.813 1.898 2.057 2.157 2.295
	300	7.883 9.358 10.552 10.482 12.435
	500	26.200 31.889 36.889 41.378 46.221
	700	65.121 75.967 90.575 104.268 116.500

Operator denotes the symmetric cartesian product defined as A B = (A × B) ∪ (B × A)


Trasanides George

Tsaousi Georgia

Adjenughwure Kingsley

Kotzampassi Katerina

Kapanidis Konstantinos

Papadopoulos Basil

Fuzzy Approach for Bibliometric Analysis of Publication Trends on Intragastric Balloon as a Minimally Invasive Procedure for Weight Loss in Obese Individuals

Keywords: Fuzzy estimators, bibliometric analysis, fuzzy sets, Bariatric Surgery, intragastric balloon, Lotka's law

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

Introduction

Obesity constitutes a worldwide epidemic socio-economic entity incurring se-rious repercussions on health status and exacerbation of obesity-related co-morbidities [START_REF] Laing | Filling the Void: A review of intragastric balloons for obesity[END_REF][START_REF] Flegal | Prevalence and trends in obesity among US adults, 1999-2008[END_REF]. These seem to exert an adverse impact on both life quality and expectancy of obese individuals, all lead-ing to spiralling health costs (Cawley J, Meyerhoefer C,2012). Although, lifestyle / behavioural and diet modification assisted by physical exercise is recommended as a first line approach to weight loss, yet numerous studies failed to support a sustainable effectiveness for the vast majority of those who attempted such practice (Norris SL et.al,2005). As to pharmacotherapy, which is considered as the second line approach, a modest efficacy associated with significant side-effects has been documented (Yanovski SZ, Yanovski JA,2014).

Thus, bariatric surgery emerges as the most reliable treatment modality for the management of severe obesity, with long-term efficacy in weight loss (Colquitt JL et al,2013; Buchwald H, Oien DM,2014). Nonetheless, the substantial cost and non-negligible surgical morbidity and mortality render bariatric surgery as a less attractive option for obese patients to pursue( Jirapinyo P, Thompson CC,2017; Buchwald H, Oien DM,2014).

To fill the void, minimally invasive non-surgical options for weight loss are gaining popularity; the intragastric balloon (IGB) being the most widely applied in clinical practice [START_REF] Laing | Filling the Void: A review of intragastric balloons for obesity[END_REF]Yorke E, et al,2016 ). The fundamental concept of an IGB is the attainment of a continuous sensation of satiety owing to mechanical gastric distension caused by an artificial gastric bezoar, thereby facil-itating decreased food intake and maintenance of a low-calorie diet (Kumar N et al 2017;Laing P et al, 2017; Yorke E, et al,2016 ). Its versatility, low invasive-ness, ease of retrieval and relative safety has resulted in the broad acceptance of its use with thousands of devices placed globally each year (Yorke E, et al,2016 In the recent literature, several bibliometric analyses have attempted to elucidate publication productivity regarding the implementation of bariatric surgery as a treatment modality of severe obesity (Aminian A, et al 2015; Ozsoy Z,& Demir E, 2017; Dabi Y eta al ,2016) or to comparably assess bariatric procedures in terms of their popularity (Ozsoy Z, &Demir E..,2018).

To the best of our knowledge, no bibliometric assessment of research output has been carried out pertaining to the literature on IGBs as a treatment modality of morbid obesity, for the 25-year period (1982 to 2017) this technique has come into clinical effect.

Thus, we conducted this bibliometric analysis on this topic using the growth of publications, publication trends over time, annual research productivity and au-thors' activity as main outcome points. Furthermore, we aimed to validate the comparable applicability of Lotka's law among authors with a single publication and those with multiple ones using a fuzzy approach.

2

Materials and Methods

Data Source and Description

The data used in this study were sourced from Scopus database. The terms Bintragastric balloons^ and Bobesity ^ with Band^ as Boolean term were used as keywords entered into the Scopus engine to retrieve data related to the objectives of this study. A search of these terms was confined to the article title, abstract and keywords of document search. The accuracy of search query was achieved with quotation marks. Documents from all subject areas, with no time period limitation and language restriction, were retrieved. Only articles published as an erratum were excluded. Search was accomplished on 20 December 2017 to avoid any bias addressed to the daily updating of the database.

Data on overall scientific production regarding the use of IGBs for obesity management were analysed for the overall growth of publications, annual re-search productivity, authors' activity, etc. To verify Lotka's law, a fuzzy approach was used.

A total of 569 articles were retrieved and further analyzed. Published articles on the use of IGBs for obesity management assessed in this study cover a time-period between 1982 up to the end of 2017. Temporal distribution of production of publications is presented in Figure 1.

Fig. 1. Timeseries of New Publications per year in the domain

Analysis of authors' productivity revealed that a total of 2001 authors contributed to the published articles. The average number of authors per article was 3.41. Detailed data on authors' activity is provided in Table 1.

Ethical approval of this study was not deemed necessary since no human sub-jects or data were involved. The first stage is the decomposition of the series and the second stage is the reconstruction of the decomposed series to get the original series. The decomposition of the series is done in two steps namely: embedding and Singular Value Decomposition (SVD). The reconstruction of the original series is done by diagonal averaging of the resulting matrix after and optional grouping of the eigenvalues step has been completed. The mathematical description of all the steps mentioned above is presented in (Golyandina, N, et al,2001;Hassani H,2007).

This technique is chosen, because it allows the timeseries to be divided into a sum of many time series with different periods and extract noise. There is no need to specify the shape of the trend beforehand.

In the second part of analysis, we try to verify, if the number publication pattern follows Lotka's law. 𝑛 𝑖 = 𝑛 1 /𝑖 𝑥 . Where 𝑛 𝑖 is the number of authors with 𝑖 publications and 𝑛 1 is the number of authors with one publication and 𝑥 is a positive real number usually assumed to be 2. In the analysis, we focus on authors with one and two publications. For this surgical procedure, there were not enough number of authors with more than 2 publications. Thus, we chose to verify 𝑛 2 = 𝑛 1 /2 𝑥 . In contrast to other studies, we do not assume that the number of publications is deterministic but we assume them to be stochastic in nature. This is a valid assumption, since we do not expect the number of publications to exactly follows Lotka's law. So, we assume that 𝑛 1 , 𝑛 2 and 𝑥 are random variables whose sampling distribution follows a normal distribution. These assumptions allow us to calculate the 95% confidence interval for 𝑛 1 , 𝑛 2 and 𝑥 using the sample estimate for the means and standard deviations of these variables. To correctly propagate the error in estimation, a fuzzy approach is adopted. In this approach, all confidence intervals, are converted to fuzzy numbers using the procedure described in (

3

Results and Discussion

Time series components

We split the time series into a sum of 15 series. This is to enable us detect any seasonality in the publications. As can be seen from the figure below, there is strong 5 years seasonality where publication number increase and decrease every 5 years. Apart from this, there are no other strong seasonality. It can be concluded that the time series is made of an increasing trend and a 5 year-periodic trend.

Fig. 2. First 6 Timeseries components as extracted by SSA technique

Trend Extraction and Fitting

The SSA technique was used to extract the trend of the original time series without specifying its form. The only parameter needed for the extraction is the window length,L which was chosen to be 15(procedure for choosing the lag is explained in the

Verification of Lotka's Law with fuzzy estimators

To verify Lotka's law, we assume that the number of authors with 1 and 2 pub-lications is stochastic. Since we do not know the true distribution, we use the sampling distribution and assume that it follows a t-distribution since the sample mean and unbiased sample standard deviation are used. We construct the 95% confidence interval for 𝑛 1 and 𝑛 2 , and transform the intervals to fuzzy numbers using the approach described in (Buckley J J,2005; Sfiris, D. S., & Papadopoulos, B. K.,2014). (Please see appendix for description of method). We then use the formula 𝑛 2 = 𝑛 1 /2 𝑥 to calculate 𝑥. Then the 95% confidence interval for 𝑥 is also calculated and transformed to a fuzzy number. The goal is to check whether for the estimated 𝑥 𝑛 1 𝑛 2 -2 𝑥 = 0

We perform the α-cut operation on the equation above as follows(see (Buckley J J,2005; Sfiris, D. S., & Papadopoulos, B. K.,2014). for arithmetic on alpha-cuts):

𝑛 ̃1 𝛼 𝑛 ̃2 𝛼 -2 𝑥 α =0
Where 𝑛 ̃1, 𝑛 ̃2,𝑥 ̃ are fuzzy numbers The resulting fuzzy numbers are shown below: The fuzzy estimator for the error in Lotka's formula shows that, publication trend indeed follows the given formula but with average x=3.0208. Using this value of x, then the formula is considered accurate to a degree of 0.95. This is a pretty high truth value which verifies that it can be used to predict number of authors with two publications. From the fuzzy number of x, it can be concluded that x=2 is highly unlikely, since its truth value equals zero. Thus, the usual Lotka's formula where x=2 does not hold in this case.

Conclusion

As the prevalence of serious obesity is constantly increased, there is a growing pressure on medical research for the evolvement of less-permanent procedures, as alternatives to surgical ones, aiming to achieve a sustainable weight loss effect that outweighs the morbidity associated with bariatric surgery. Use of IGBs emerges as a technically easier and of minimal-risk procedure en-suring a significant weight loss in selected obese patients (Kotzampassi K et al,2012) . IGBs being currently in practice are delivered endoscopically into the stomach to reduce gastric capacity and inflated with either saline or air. It has been advocated as an effective method of inducing adequate weight loss in non-morbid obesity for whatever reason or aesthetic purposes, prior to a more definitive treatment through bariatric surgery in an attempt to minimize complication rates and technical difficulties of surgery or as an adjunct to conservative treatment in the morbidly obese denied anaesthesia and/or surgery (Fernandes M et al,2007, Kotzampassi K et al,2006, Kotzampassi K et al,2008). Typically, IGB devices are removed after 6 months, as the risk of balloon deflation and possible migration is heightened substantially beyond this time-period (Papavramidis TS et al,2012). Bibliometric analysis is regarded as one scientometry tool used to assess the scientific activities, according to the four major domains of published articles, such as authors, citations, references and publications. In the present study we decided to evaluate both the trend of publications and alos verify Lotka's law using a fuzzy approach. The analysis shows that publication trend follows Price's law. The extracted trend is exponential and in accordance with what is found in other studies. Lotka's formula using one and two publications, is found to hold but for mean x=3.025 . Thus using the original value x=2 will grossly over-estimate number of authors with two publications. Further analysis is needed to verify, if this hold for higher number of publications also. The proposed approach also has some limitations regarding the calculations with alphacuts. Since the exact membership functions are not known, the me membership values are approximated.

  Buckley J J,2005; Sfiris, D. S., & Papadopoulos, B. K.,2014) A description of the fuzzy theory and fuzzy estimators is presented in (Buckley J J,2005; Sfiris, D. S., & Papadopoulos, B. K.,2014).

  (Golyandina, N, et al,2001; Hassani H,2007). The extracted trend shows an exponential increase in number of publications as expected following Price's law. An exponential function is fitted to the observed data. The fitted function is in line with Prices' law and corresponds to those results found in(Dabi Y, et al,2016).

Fig. 3 .

 3 Fig. 3. Price's Law exponential fitting of publication trend

Fig. 4 .Fig. 5 . 1 Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 4516789 Fig. 4. Fuzzy estimation for number of authors with one publication (mean value has a truth value of 1)

Fig. 10 .

 10 Fig. 10. Publication trend for authors with 2 publications showing observed and estimated trend for x=2, and x=3.0208

  ). Consequently, numerous publications have assessed the applicability of various types of IGBs from their inception to clinical practice up to nowadays(Laing P et al, 2017; Yorke E, et al,2016; Kotzampassi K et al,2012; Kotzampassi K et al,2014) . Bibliometric analysis incorporates statistical and mathematical methods with a view to provide a quantitative and qualitative analysis of scientific evidence, as well as the course of research performance and productivity on a certain scien-tific topic and period of time(Ellegaard O, Wallin JA ,2015; Mishra L et. al,2016;,Thompson DF, Walker CK,2015). Bibliometric indicators have been extensively used in various medical disciplines to outline research trends, sug-gest future research ideas and challenge public health issues(Khan NR et al,2014;. Sweileh WM et al ,2015, Zyoud SH et al,2015 )

Table 1 .

 1 Authors' Activities

		No of authors	No of articles	Percentage %
		1533	1	76,7
		466	> 1	23,3
		290	2	14,5
		73	3	3,65
		95	4 to 9	4,75
		8	> 10	0,4
			11 to 27
	2.2	Data Analysis Techniques

This analysis has two parts, first we look at the publication trend by analysing the time series of publications. The trend is extracted from the original series using the singular spectrum analysis technique

(Golyandina, N, et al,2001; Has-sani H,2007)

. Singular Spectrum Analysis (SSA) is a relatively new technique for time series decomposition, reconstruction and prediction. The SSA is a two-stage process.

Appendix

Let 𝑋 be a universal set. A fuzzy set is defined as any function of the form 𝐴: 𝑋 → [0,1]. 𝐴(𝑥) is the membership degree of 𝑥 in the fuzzy set 𝐴. The 𝛼-cuts of a fuzzy set 𝐴 are defined as the sets:

The 𝛼-cuts of a fuzzy number 𝐴 can written as intervals of the form 𝐴 = [𝐴

Non-Asymptotic fuzzy estimators

Buckely 2005 introduced fuzzy estimators. Papadopoulos and Sfiris,2014 generalize fuzzy fuzzy estimators and showed new ways to construct them.

Let 𝑋 1 , 𝑋 2 , … 𝑋 𝑛 be a random sample of size 𝑛 from a distribution with unknown parameter 𝜃 and let [𝜃 1 (𝛼) , 𝜃 2 (𝛼)] denote the (1 -𝛼)100% confidence intervals for 𝜃. Using any monotonic, continuous and onto function:

, 0.5],𝛾 ∈ (0,1), then a family of fuzzy estimators is derived by the following super-imposed intervals: 𝜃 α 𝛾 =[𝜃 1 (2ℎ(𝛼)) , 𝜃 2 (2ℎ(𝛼))], 𝛼 ∈ (0,1] It is shown that these fuzzy estimators are non-asymptotic and proper fuzzy numbers with compact support. Note that the explicit membership function can also be derived as :

𝜃 ̃𝛾(𝑢) = 𝑚𝑖𝑛{𝜃 1 -1 (𝑢) , -𝜃 2 -1 (-𝑢), 1 }

For the confidence interval of the mean, this membership function can be derived as:

The α-cut of this fuzzy number is :