David Griol
email: dgriol@inf.uc3m.es

José Antonio Iglesias
email: jiglesia@inf.uc3m.es

Agapito Ledezma
email: ledezma@inf.uc3m.es

Araceli Sanchis

A practical application of evolving fuzzy-rule-based classifiers for the development of spoken dialog systems ⋆

Keywords: Spoken Dialog Systems, Evolving Classifiers, Spoken Human-Machine Interaction, Statistical Methodologies, Dialog Management

In this paper we present a statistical approach based on evolving Fuzzy-rule-based (FRB) classifiers for the development of dialog managers for spoken dialog systems. The dialog managers developed by means of our proposal select the next system action by considering a set of dynamic fuzzy rules that are automatically obtained by means of the application of the FRB classification process. Our approach has the main advantage of taking into account the data supplied by the user throughout the complete dialog history without causing scalability problems, also considering confidence measures provided by the recognition and understanding modules. The use of EFS also allows to process streaming data on-line in real time, thus dynamically evolving the structure and operation of the dialog model based on the interaction of the dialog system with its users. We also describe the application of our proposal to develop a dialog system providing railway information.

Introduction

Spoken Dialog Systems (SDSs) are computer programs that receive speech as input and generate synthesized speech as output, engaging the user in a dialog that aims to be similar to that between humans [START_REF] Pieraccini | The Voice in the Machine: Building Computers Understand Speech[END_REF]. These interfaces make technologies more usable, as they ease interaction, allow integration in different environments, and make technologies more accessible, especially for disabled and elderly people [START_REF] Vipperla | Spoken dialogue interfaces for older people[END_REF].

Usually, SDSs carry out five main tasks [START_REF] Pieraccini | The Voice in the Machine: Building Computers Understand Speech[END_REF][START_REF] Schatzmann | A Survey of Statistical User Simulation Techniques for Reinforcement-Learning of Dialogue Management Strategies[END_REF][START_REF] Griol | A Statistical Approach to Spoken Dialog Systems Design and Evaluation[END_REF]: Automatic Speech Recognition (ASR), Spoken Language Understanding (SLU), Dialog Management (DM), Natural Language Generation (NLG), and Text-To-Speech Synthesis (TTS). These tasks are typically implemented in different modules of the system's architecture.

Although dialog management is only a part of the development cycle of spoken dialog systems, it can be considered one of the most demanding tasks given that it decides the next action of the system and encapsulates the logic of the speech application [START_REF] Wilks | Some background on dialogue management and conversational speech for dialogue systems[END_REF]. The selection of a specific system action depends on multiple factors, such as the output of the speech recognizer (e.g., measures that define the reliability of the recognized information), the dialog interaction and previous dialog history (e.g., the number of repairs carried out so far), the application domain (e.g., guidelines for customer service), knowledge about the users, and the responses and status of external back-ends, devices, and data repositories. Given that the actions of the system directly impact users, the dialog manager is largely responsible for user satisfaction. This way, the design of an appropriate dialog management strategy is at the core of dialog system engineering.

As an attempt to reduce the time and effort required for system implementation and carry out rapid system prototyping, statistical approaches for dialog management are gaining increasing interest. These approaches enable automatic learning of dialog strategies, thus avoiding the time-consuming process that hand-crafted dialog design involves. Statistical models can be trained from real dialogs, modeling the variability in user behaviors. Although the construction and parameterization of these models depend on expert knowledge about the task to be carried out by the dialog system, the final objective is to develop systems that are more robust for real-world conditions, and that are easier to adapt to different users and tasks [START_REF] Schatzmann | A Survey of Statistical User Simulation Techniques for Reinforcement-Learning of Dialogue Management Strategies[END_REF].

In this paper we present a statistical approach for the development of dialog managers, in which the next system prompt is selected by means of a novel approach for online classifying based on Evolving Fuzzy Systems (EFS) [START_REF] Angelov | Evolving fuzzy-rule-based classifiers from data streams[END_REF]. These systems have an evolving structure which is updated according to the input samples. A prototype for this kind of classifiers is a data sample that groups several samples. The classifier is initialized with the first data sample and then, each data sample is classified into one of the existing prototypes. Finally, based on the potential of the new data sample to become a prototype, form a new prototype, or replace an existing one. Unlike other statistical approaches, our approach has the advantage of taking into account the data supplied by the user throughout the complete dialog without causing scalability problems. Confidence measures provided by the recognition and the understanding modules are also taken into account in the selection of the next system response.

The most widespread methodology for machine-learning of dialog strategies consists of modeling human-computer interaction as an optimization problem using Markov Decision Processes (MDP) and reinforcement methods [START_REF] Levin | A stochastic model of human-machine interaction for learning dialog strategies[END_REF]. The main drawback of this approach is that the large state space of practical spoken dialog systems makes its direct representation intractable. As described in [START_REF] Young | The Hidden Information State Approach to Dialogue Management[END_REF], Partially Observable MDPs (POMDPs) outperform MDP-based dialog strategies since they provide an explicit representation of uncertainty. This enables the dialog manager to avoid and recover from recognition errors by sharing and shifting probability mass between multiple hypotheses of the current dialog state. Other interesting approaches for statistical dialog management are based on modeling the system by means of stochastic Finite-State Transducers [START_REF] Planells | An Online Generated Transducer to Increase Dialog Manager Coverage[END_REF], or Bayesian Networks [START_REF] Meng | The Use of Belief Networks for Mixed-Initiative Dialog Modeling[END_REF].

The use of EFS allows us to cope with huge amounts of data, process streaming data on-line in real time, and evolve the structure of a dialog model that defines an activity based on the human-computer interaction. The dialog model is designed and treated as a changing model which constantly reflects the changes in the way the dialog system interacts with the users. Thus, it can start learning "from scratch" since the number of fuzzy rules do not need to be defined a priori. However, EFS can also be used to evolve an initial dialog model which has been previously defined taking into account expert knowledge about the specific interaction domain. In addition, the proposed solution addresses other important related requirements such as processing huge amounts of sequential data in an optimized way.

After this introduction, the remainder of the paper is organized as follows. Section 2 describes the proposed methodology for dialog management and its specific application for the eclass0 classifier. This section also presents a practical application of our proposal to develop a spoken dialog system providing railway information and the results of a preliminary evaluation. Section 3 presents the conclusions and suggests some future work guidelines.

Our proposed methodology for dialog management

The process followed by a statistical dialog manager for the selection of the next system response is summarized in Figure 1 [START_REF] Young | The Statistical Approach to the Design of Spoken Dialogue Systems[END_REF]. As this figure shows, the user has an internal state S u corresponding to a goal that is trying to accomplish and the dialog state S d represents the previous history of the dialog. Based on the user's goal prior to each turn, the user decides some communicative action (also called intention) A u , expressed in terms of dialog acts and corresponding to an audio signal Y u . Then, the speech recognition and language understanding modules take the audio signal Y u and generate the pair ( Ãu , C).

This pair consists of an estimate of the user's action A u and a confidence score that provides an indication of the reliability of the recognition and semantic interpretation results. This pair is then passed to the dialog model, which is in an internal state S m an decides what action A m the dialog system should take. This action is also passed back to the dialog manager so that S m may track both user and machine actions. The language generator and the text-to-speech synthesizer take A m and generate an audio response Y m . The user listens to Y m and attempts to recover A m . As a result of this process, users update their goal state S u and their interpretation of the dialog history S d . These steps are then repeated until the end of the dialog.

In our statistical approach for dialog management, we propose that, given a new user turn, the statistical dialog model makes the assignation of a system response A m according to the result of a classification process. The classification function can be defined in several ways. We propose the use of the eClass (evolving Classifier) family. During the training of these classifiers, a set of fuzzy rules Fig. 1. Set of internal processes in a spoken dialog system that describes the most important observed features for the classification of each class (i.e., system prompt) is formed. These rules can also be constantly adjusted to the available training data. One of the advantages of eClass is that it does not require parameter optimization as its only parameter 'scale' can be directly inferred from the training data. This technique [START_REF] Angelov | An approach to online identification of Takagi-Sugeno fuzzy models[END_REF] is based on partitioning the data space into overlapping local regions through Recursive Density Estimation (RDE) and associating clusters (respectively fuzzy sets) to them.

The eClass family includes two different architectures and on-line learning methods:

-eClass0 with the classifier consequents representing class label [START_REF] Angelov | An approach to online identification of Takagi-Sugeno fuzzy models[END_REF].

-eClass1 for regression over the features using first order eTS fuzzy classifier.

In this paper, we focus on the eClass0 classifier. eClass0 possesses a zeroorder Takagi-Sugeno consequent, so a fuzzy rule in the eClass0 model has the following structure:

Rule i = IF (F eature 1 is P 1 ) AN D . . . . . . AN D (F eature n is P n ) T HEN Class = c i (1)
where i represents the number of rule; n is the number of input features (observations corresponding to the different attributes and concepts defined for the semantic representation of the user's utterances); the vector F eature stores the observed features, and the vector P stores the values of the features of one of the prototypes of the corresponding class c i ∈ {set of different classes}. Each class is then associated to a specific system action (response).

The dialog manager must consider the concepts and values for the attributes provided by the user throughout the previous history of the dialog to select the next system action. For the dialog manager to take this decision, we have assumed that the exact values of the attributes are not significant. They are important for accessing databases and for constructing the output sentences of the system. However, the only information necessary to predict the next action by the system is the presence or absence of concepts and attributes. Therefore, the codification we use for each concept and attribute provided by the SLU module is in terms of three values, {0, 1, 2}, according to the following criteria:

-(0): The concept is unknown or the value of the attribute is not given; -(1): The concept or attribute is known with a confidence score that is higher than a given threshold; -(2): The concept or attribute has a confidence score that is lower than the given threshold.

The eClass0 model is composed of several fuzzy rules per class (the number of rules depends on the heterogeneity of the input data of the same class). During the training process, a set of rules is formed from scratch using an evolving clustering approach to decide when to create new rules. The inference in eClass0 is produced using the "winner takes all" rule and the membership functions that describe the degree of association with a specific prototype are of Gaussian form.

The potential (Cauchy function of the sum of distances between a certain data sample and all other data samples in the feature space) is used in the partitioning algorithm. However, in these classifiers, the potential (P) is calculated recursively (which makes the algorithm faster and more efficient). The potential of the k th data sample (x k ) is calculated by means of equation 2 [START_REF] Angelov | Evolving fuzzy-rule-based classifiers from data streams[END_REF]. The result of this function represents the density of the data that surrounds a certain data sample.

P (x k ) = 1 1 + ∑ k-1 i=1 distance(x k ,xi) k-1 ( 2
)
where distance represents the distance between two samples in the data space. The potential can be calculated using the euclidean or the cosine distance. In this case, cosine distance (cosDist) is used to measure the similarity between two samples; as it is described in equation 3.

cosDist(x k , x p ) = 1 - ∑ n j=1 x kj x pj √ ∑ n j=1 x 2 kj ∑ n j=1 x 2 pj ( 3
)
where x k and x p represent the two samples to measure its distance and n represents the number of different attributes in both samples.

Note that the resolution of equation 2 requires all the accumulated data sample available to be calculated, which contradicts to the requirement for realtime and on-line application needed in the proposed problem. For this reason, in [START_REF] Angelov | Evolving fuzzy-rule-based classifiers from data streams[END_REF] it is developed a recursive expression for the cosine distance. The proposed formula is as follows:

P k (z k ) = 1 2 - 1 (k-1) √ ∑ n j=1 (z j k ) 2 B k ; k = 2, 3... where : B k = n ∑ j=1 z j k b j k ; b j k = b j (k-1) + √ (z j k ) 2 ∑ n l=1 (z l k ) 2 and b j 1 = √ (z j 1 ) 2 ∑ n l=1 (z l 1 ) 2 ; j = [1, n + 1]; P 1 (z 1 ) = 1 (4)
where z k represents the k th data sample (x k ) and its corresponding label (z = [x, Label]). Using this expression, it is only necessary to calculate (n+1) values where n is the number of different subsequences obtained; this value is represented by b, where b j k , j = [1, n] represents the accumulated value for the k th data sample.

In this case, a specific system action can be represented by several rules, depending on the heterogeneity of the samples that represent the same action. Thus, a class could be represented by one or several prototypes. The different prototypes that represent a system action are obtained from the input data and they are updated constantly. However, an initial rule-based model can be defined (if necessary) by hand as start point of the classifier. In this sense, new prototypes are created or existing prototypes are removed if necessary.

Practical application: a spoken dialog system providing railway information

We have applied our proposal to develop a mixed-initiative spoken dialog system to provide railway information system using spontaneous speech in Spanish. The system integrates the CMU Sphinx-II system speech recognition module1 . As in many other conversational agents, the semantic representation chosen for dialog acts of the SLU module is based on the concept of frame [START_REF] Minsky | A Framework for Representing Knowledge[END_REF]. This way, one or more concepts represent the intention of the utterance, and a sequence of attribute-value pairs contains the information about the values given by the user.

For the task, we defined eight concepts and ten attributes. The eight concepts are divided into two groups:

1. Task-dependent concepts: they represent the concepts the user can ask for (Timetables, Fares, Train-Type, Trip-Time, and Services). 2. Task-independent concepts: they represent typical interactions in a dialog (Acceptance, Rejection, and Not-Understood).

The attributes are: Origin, Destination, Departure-Date, Arrival-Date, Class, Departure-Hour, Arrival-Hour, Train-Type, Order-Number, and Services.

A total of 51 system responses were defined for the task (classified into confirmations of concepts and attributes, questions to require data from the user, and answers obtained after a query to the database).

Using the previously described codification for the concepts and attributes, when a dialog starts (in the greeting turn) all the values are initialized to "0". The information provided by the users in each dialog turn is employed to update the previous values and obtain the current ones, as Figure 2 shows.

This figure shows the semantic interpretation and confidence scores (in brackets) for a user's utterance provided by the SLU module. In this case, the confidence score assigned to the attribute Date is very low. Thus, a "2" value is added in the corresponding position for this attribute. The concept (Hour ) and the attribute Destination are recognized with a high confidence score, adding a "1" value in the corresponding positions. The set of features for the classifier includes the codification of the different concepts and attributes that can be provided by the user and the taskindependent information provided in the last user turn (none in this case). A total of 49 rules for the task were obtained with eClass0. Figure 3 shows the structure of these rules. Using them, the dialog manager would select the class ′ SystemResponse23 ′ , which corresponds to a system confirmation of the departure date. This process is repeated to predict the next system response after each user turn.

An initial corpus of 900 dialogs (10.8 hours) was acquired for the task by means of the Wizard of Oz technique with 225 real users, for which an initial dialog strategy was defined by experts [START_REF] Griol | A Statistical Approach to Spoken Dialog Systems Design and Evaluation[END_REF]. A set of 20 scenarios was used to carry out the acquisition. Each scenario defined one or two objectives to be completed by the user and the set of attributes that they must provide. The corpus consists of 6,280 user turns, with an average number of 7.7 words per turn. The corpus was split into a training subset of 4,928 samples (80% of the corpus) and a test subset of 1,232 samples (20% of the corpus).

We defined the following three measures to evaluate the statistical dialog model: i) Matching: percentage of responses provided by the system that match exactly the initial dialog strategy defined by the experts; ii) Coherence: percentage of responses that are coherent with the current dialog state but do not match the initial strategy; and iii) Error : percentage of responses that could cause a dialog failure.

We have test the behavior of our proposal comparing it with different definitions of the classification function used to determine the next system response. In this work, we have used three approaches for the definition of the classification function: a multilayer perceptron (MLP), a multinomial naive Bayes classifier, F RB -RailwayT ask(eClass0) : 3. Set of rules for the dialog manager obtained with eClass0 for the railway task and finite-state classifiers. We also defined three types of finite-state classifiers: bigram models, trigram models, and Morphic Generator Grammatical Inference (MGGI) models [START_REF] Segarra | Construction of Language Models using Morfic Generator Grammatical Inference MGGI Methodology[END_REF].

IF (T imetables is 1) AN D (F ares is 0) AN D • • • AN D (N ot -U nderstood is 0) T HEN Class = ′ Ask -Date ′ IF (T imetables is 2) AN D (F ares is 0) AN D • • • AN D (N ot -U nderstood is 1) T HEN Class = ′ Conf irm -T imetables ′ IF (T imetables is 0) AN D (F ares is 1) AN D • • • AN D (N ot -U nderstood is 0) T HEN Class = ′ P rovide -F ares ′ • • • IF (T imetables is 1) AN D (F ares is 2) AN D • • • AN D (N ot -U nderstood is 1) T HEN Class = ′ Close -Dialog ′ Fig.
Table 1 shows the results obtained. As it can be observed, the Fuzzy-rulebased classifier provides satisfactory results in terms of the percentage of correct responses selected (Matching and Coherence measures) and responses that could cause the failure of the dialog (Error measure). With regard the rest of classifiers, the MLP classifier is the one providing the closest results to our proposal. The table also shows that among the finite-state model classifiers, the bigram and trigram classifiers are worse than the MGGI classifier, this is because they cannot capture long-term dependencies. The renaming function defined for the MGGI classifier seems to generate a model with too many states for the size of the training corpus, therefore, this classifier could be underestimated.

Dialog manager

Matching Secondly, we have evaluated our proposal with the acquisition of 100 dialogs by means of a user simulator [START_REF] Griol | A Statistical Approach to Spoken Dialog Systems Design and Evaluation[END_REF]. We considered the following measures: i) Dialog success rate (Success); ii) Average number of turns per dialog (nT ); iii) Confirmation rate (Confirmation); and iv) Error correction rate (ECR). The confirmation rate was computed as the ratio between the number of explicit confirmation turns and the total number of turns in the dialog. The ECR was computed as the number of errors detected and corrected by the dialog manager divided by the total number of errors.

The results presented in Table 2 show that in most cases the automatically learned dialog model has the capability of correctly interacting with the user. The dialog success depends on whether the system provides the correct data for every objective user's query. All of the objectives defined are achieved in 93.5% of the dialogs. The analysis of the main problem detected in the acquired dialogs shows that, in some cases, the system did not detect the introduction of data with a high confidence value due to errors generated by the ASR that were not detected by the dialog manager. However, the evaluation confirms a good operation of the approach since the information is correctly given to the user in the majority of cases. The confirmation and error correction rates have also a remarkable impact on the described system performance.

Conclusions and Future Work

In this paper, we have presented a statistical methodology for the development of dialog managers and the optimization of dialog strategies in spoken dialog systems. The selection of the following system response is based on a classification process that takes into account the history of the dialog. The most important contribution of our work consists of the use of Evolving Fuzzy Systems (EFS) to complete this classification. As a result of the application of our proposal, the dialog model is modeled by a set of automatically obtained dynamic rules that are applied to select the next system response. The use of EFS allows us to cope with huge amounts of data, and process streaming data on-line in real time.

We have described a practical application of our proposal to develop and evaluate a spoken dialog system providing railway information. A codification of the information sources has been proposed to facilitate the correct operation of the eClass0 classification function. This representation allows the system to automatically generate a specialized answer that takes the current situation of the dialog into account. Task-dependent information is isolated from the model taking into account whether the user has provided a given piece of information related to the task and also the confidence scores assigned by the ASR and NLU modules.This allows not only to cope with the situations observed the training corpus, but also to manage unseen situations by selecting the most convenient system action. The results of the evaluation shows the correct operation of the learned dialog manager with regard other definitions of the classification functions.

Future work will be oriented to deploy and evaluate our proposal in additional domains. As the dialog model is learned from a corpus of training samples, the performance of the dialog manager depends on the quality and size of the corpus used to learn the model. For this reason, we also want to evaluate the influence of the main features of the training corpus in the quality of the dialog model obtained by means of our proposal.

Fig. 2 .

 2 Fig. 2. Excerpt of a dialog with its correspondent representation of the task-dependent and active task-independent information for one of the dialog turns

Table 1 .

 1 Results of the evaluation of the different classification functions

	Coherence Error

Table 2 .

 2 Results of the objective evaluation with real users

	Success nT Confirmation ECR

cmusphinx.sourceforge.net

⋆ This work has been supported by the Spanish Government under i-Support (Intelligent Agent Based Driver Decision Support) Project (TRA2011-29454-C03-03).


Josef Hynek
email: josef.hynek@uhk.cz

Sequence Matching Genetic Algorithm for Square Jigsaw Puzzles

Keywords: genetic algorithms, image reconstruction, pattern matching, puzzle solving

Our paper presents a new method for solving the rectangle piece jigsaw puzzle problem. The puzzle image is RGB full color and because of uniform shape of the individual pieces the process of puzzle assembly is based on information of the pixel values along the border line of the piece only. We have utilized a genetic algorithm that searches for the optimal piece arrangement using dissimilarity between adjacent pieces as the measure of progress. Unlike the previous attempts to utilize genetic algorithms to solve the problem, we have proposed a new heuristic asexual operator that aims at identification of points of fraction within partially assembled picture, extraction of supposed sequence of correctly joint pieces, and its insertion into a new position in such a way that, if possible, the segment is enlarged. Our approach has been successfully tested and the algorithm is capable of solving puzzles consisting of several hundred pieces.

INTRODUCTION

Solving puzzles is usually viewed as a special example of pattern matching recognition problem. At the first sight it might be viewed as a pleasant and very popular widespread game for players of all ages but we should not forget the complexity of jigsaw puzzles. It is a well-known fact that jigsaw puzzles and their equivalents belong to NP-complete class [START_REF] Demaine | Jigsaw Puzzles, Edge Matching, and Polyomino Packing: Connections and Complexity[END_REF] and thus provide serious obstacles for various attempts to solve them automatically. Moreover, there are many practical applications ranging from image mosaicking up to the reconstruction of archaeological artifacts and assembly of shredded documents [START_REF] Richter | Learning to Reassemble Shredded Documents[END_REF]. The first jigsaw puzzle was produced in 1760 in London [START_REF] Gindre | A Criterion-based Genetic Algorithm Solution to the Jigsaw Puzzle NP-Complete Problem[END_REF] and then the game of puzzle gradually expanded all over the world. There are various types of jigsaw puzzle and an integrated classification terminology is missing. However, there are some basic classification terms that are commonly used and we will mention here at least the most important of them. The term standard jigsaw puzzle is usually used to describe puzzles created by cutting pictures printed on firm stiff paper into interlocking patterns of pieces. The pieces have a distinctive geometrical shape that provides an important piece of information that is used together with image information in the process of puzzle assembly. It is necessary to point out that there are also so cold apictorial types of puzzle that do not have any picture not chromatic information and the shape of the pieces is the only clue to solve the problem.

We are interested in square jigsaw puzzles where the individual pieces lack characteristic curvilinear shape because all the pieces have straight borders and uniform rectangular size. In this case pictorial information provides us only information and guidance that is available there. This type of puzzle is sometimes called as edgematching puzzle [START_REF] Demaine | Jigsaw Puzzles, Edge Matching, and Polyomino Packing: Connections and Complexity[END_REF] in order to emphasize that the adjacent pieces match along their respective edges. These puzzles are very challenging as there is no guarantee that two pieces that fits together should be together and the only way to proof the correctness of individual matches is the completion of the entire solution.

Of course, the simple categorization of different kinds of puzzle above is very superficial and definitely incomplete, but for the sake of this paper we do not need to go further in this direction. For more details concerning different kinds of puzzle and its classification please refer to [START_REF] Gindre | A Criterion-based Genetic Algorithm Solution to the Jigsaw Puzzle NP-Complete Problem[END_REF].

The rest of the paper is organized as follows. Section 2 summarizes various approaches to the automatic solution of jigsaw puzzles with the particular focus on utilization of genetic algorithms. Section 3 is the key part of our paper and the newly proposed solution is described there. Our experimental results are briefly described in Section 4 and it is followed by discussion and conclusions.

Previous Work

Many attempts to utilize various problem solving techniques were made since 1964 when the first paper describing computer solution of apictorial jigsaw puzzle was published [START_REF] Freeman | Apictorial Jigsaw Puzzles: The Computer Solution of a Problem in Pattern Recognition[END_REF]. Is has been pointed out that the difficulty encountered when attempting to program a computer to solve jigsaw puzzles relates to three different aspects of the problem: 1) the description of the pieces, 2) the manipulation of the pieces (rotating and matching), and 3) the evaluation of the correct matching of the individual pieces. The complexity of the problem implies that a brute-force approach, which may work for small puzzles, becomes impractical as the number of pieces is increased [START_REF] Freeman | Apictorial Jigsaw Puzzles: The Computer Solution of a Problem in Pattern Recognition[END_REF].

It is necessary to emphasize that majority of attempts that dates more that fifteen years back focused on apictorial puzzles. It means that the relevant experiments are concerned with geometric shape information of the puzzle pieces and chromatic (pictorial) information is neglected (see, for example, [START_REF] Wolfson | Solving jigsaw puzzles by computer[END_REF] or [START_REF] Webster | Isthmus critical points for solving jigsaw puzzles in computer vision[END_REF]). These articles brought up gradually improved methods and techniques for representation of boundaries of puzzle pieces and their subsequent efficient pertaining and correct matching.

Perhaps the most efficient algorithm for automatic solution of jigsaw puzzles was presented by Goldber et al. [START_REF] Goldberg | A global approach to automatic solution of jigsaw puzzles[END_REF]. It is based on the similar ideas as the already mentioned algorithm described in [START_REF] Wolfson | Solving jigsaw puzzles by computer[END_REF], but it utilizes better and more global matching techniques. In particular, use of fiducial points to align adjacent pieces of puzzle together with optimization of partial solutions based on the method of global relaxation worked surprisingly well and allowed authors to solve the puzzles consisting of more than two hundred pieces [START_REF] Goldberg | A global approach to automatic solution of jigsaw puzzles[END_REF].

An obvious way to improve efficiency and robustness of the algorithm is a suitable utilization of additional pieces of information on the top of a distinctive geometric shape of the puzzle pieces. Colours or textures on the boundaries of these pieces could provide the algorithm further clues and direction if such information is available (see [START_REF] Chung | Jigsaw puzzle solver using shape and color[END_REF] or [START_REF] Makridis | An innovative algorithm for solving jigsaw puzzles using geometrical and color features[END_REF]). However, as we have mentioned above, there are also so called edgematching puzzles where all the pieces have straight borders and uniform rectangular size and pictorial information provides us the only information and guidance that can be used in order to solve it.

Toyama et al. [START_REF] Toyama | Assembly of puzzles using a genetic algorithm[END_REF] proposed a method for solving the rectangle piece jigsaw puzzle problem using a genetic algorithm. The picture of their puzzle is painted only in black and white and so the pieces of the puzzle are represented as binary images. The assembly of the puzzle is then performed using information of the pixel value on the border line of the adjacent pieces. They utilised the populations of 200 randomly generated individuals representing candidate solutions and employed two genetic operators (2-point crossover and self-crossover) to produce offspring and run the evolutionary process. They reported that the method described there correctly assembled all pieces in the 8x8-piece puzzle.

This team of researchers collaborated on another interesting paper [START_REF] Murakami | Assembly of puzzles by connecting between blocks[END_REF] and developed an improved method for solving the rectangle piece jigsaw puzzle assembly problem, but no evolutionary algorithm is used in this process any more. The assembly of the puzzle is once again performed only using information of the pixel value on the border line of the pieces and this time a puzzle image is RGB full colour. Pieces are connected by a matching function between two pieces and a simple method connects together a single piece to a block that is defined as a group of already connected pieces. According their results the proposed method correctly assembled all pieces in 16x12-piece puzzles.

Alajlan [START_REF] Alajlan | Solving Square Jigsaw Puzzles Using Dynamic Programming and the Hungarian Procedure[END_REF] experimented with gray picture puzzles and he used dynamic programming to facilitate non-rigid alignment of border pixels for local matching of the puzzle pieces. Moreover, instead of the classical best-first search, his algorithm simultaneously positioned the neighbours of a puzzle pieces during the search using the socalled Hungarian procedure. This procedure begins with a starting piece of the puzzle and then locates four adjacent pieces of puzzle to the initial piece in the way the sum of border distances is minimal. In order to make the algorithm robust, every puzzle piece was considered as starting piece at various starting locations. Experiments using several images demonstrated that the proposed algorithm correctly assembled puzzles up to 8×8 pieces. However, the author of this article acknowledged that its performance deteriorates as the number of pieces exceeds 64.

Finally, we would like to mention one completely different approach that was described by Gindre et al. [START_REF] Gindre | A Criterion-based Genetic Algorithm Solution to the Jigsaw Puzzle NP-Complete Problem[END_REF]. Their algorithm is a part of the development of Intelligent Robotic System that solves an unknown jigsaw puzzle. The system uses pattern recognition techniques as edge and feature detection in conjunction with genetic algorithms. The novelty of this as approach is based on a completely different encoding scheme. The candidate solutions are represented by relevant graphs describing the interconnections between puzzle pieces and that is why the corresponding chromosomes are defined by the relevant adjacency matrix.

The common limitation of the above presented approaches is that just small scale puzzles could be tackled and solved. That is the fact that fostered our interest in this area and based on a detailed analysis of the methods that were already tested we have designed a hybrid genetic algorithm for jigsaw puzzle problem.

Problem Definition and Representation

We have shown that there are various types of jigsaw puzzles and therefore it is necessary to start with its definition. We have decided to experiment with rectangular pictorial puzzle as it is described in [START_REF] Toyama | Assembly of puzzles using a genetic algorithm[END_REF]. Therefore, we assume that it consists of NxM puzzle pieces that do not rotate and the assembly of the puzzle is performed by using the pixel values on the border lines of the individual pieces. Whereas Toyama worked with black and white pictures only where the pixel values are 0 and 1 respectively, we have decided to use color 24-bit RGB images. Despite the fact that the comparison of border lines is very easy and straightforward, there is no guarantee that there is an absolute correspondence between border lines of the neighboring pieces. It depends on the particular picture and irregularities along the border make the problem difficult. It is clear that algorithm could be even mislead in the situation when wrong piece of the puzzle is regarded as more suitable from the point of borders matching than the right one. Hence, it is evident that the problem cannot be solved by the simple local piece matching. The example of the puzzle used in our experiments is shown on Fig. 1.

Having defined the problem we can decide what kind of representation of partial (or candidate) solutions would be suitable for our algorithm. We have realized that graph representation used by Gindre et al. [START_REF] Gindre | A Criterion-based Genetic Algorithm Solution to the Jigsaw Puzzle NP-Complete Problem[END_REF] is rather inconvenient for genetic operators and the pertinent manipulations and that is why we have chosen simple NxM matrix for each individual and each cell c i,j corresponds to one piece of puzzle. The matrix depicts the relevant arrangement of puzzle pieces in a very straightforward manner that the left/right/up/down adjacency in the matrix means the same interrelationship between the relevant puzzle pieces. This representation is very convenient for fitness evaluation as well as for utilization of genetic operators.

Each piece of the puzzle is represented by matrix LxLx3, where L is the width and length of the piece (i.e. the number pixels along the border). The fitness function of our algorithm is based on two measures of dissimilarity. Horizontal dissimilarity between two pieces x i and x j , where x j is placed to the right of x i is defined as

, )) , 1 , ( ) , , ( ( ) , ( 1 3 1
2      L l k j i j i h c l x c L l x x x D (1)
where c stands for intensity of red, green, and blue color and each value is integer from 0 to 256. Likewise, vertical dissimilarity between two pieces x i and x j , where x i is placed on the top of x j is defined as

, )) , , 1 ( ) , , ( ( ) , ( 1 3 1 2      L l k j i j i v c l x c l L x x x D (2)
where square difference between pixels of the last row of x i and the first row of x j is computed. These measures of dissimilarity will be calculated many times during the process of evaluation of thousands of individuals and therefore it is advantageous to create a lookup table covering all the relevant values at the beginning of run. The fitness function of each individual is calculated as the reverse value of the sums of horizontal and vertical differences of the relevant puzzle pieces. Providing that individual a is represented by NxM matrix C as we have described it above, then its fitness is defined as

. ) , ( ) , ( 1
1 ) ( 1 1 1 , 1 , 1 1 1 1 , ,               N i M j j i j i v N i M j j i j i h c c D c c D a f ( 3
)
It is evident that the important feature of fitness function f(a) is that individual with smaller sum of dissimilarities is fitter and thus will have a higher chance of being selected. We have employed standard biased roulette wheel selection method together with an elitist approach when the best individuals (10%) from current population are copied without any change to the newly created population.

The novelty of our approach is based on the utilization of problem specific operator. Besides the common partially-matched crossover (PMX) that assures legitimacy of created offspring we have employed a new operator that helps to preserve the already created segments of the picture being assembled. The main drawback of PMX is that crossover points are selected randomly and therefore it is easy to split a sequence of correctly assembled pieces. Nevertheless, it provides us necessary variability, exchanges pieces of information between different individuals and thus it is important for the whole process of evolution.

We have experimented with asexual heuristic operator that is designed to enlarge the already created segments of the picture. These segments are not only scattered over the whole candidate solution but they are usually located at incorrect locations. We are unable to assign them to their correct place within the picture as the algorithm does not use this piece of information, nevertheless we can try to connect them with other segments believing that larger segments will finally settle at the right positions.

First of all, we identify likely points of fraction where the horizontal measure of dissimilarity is relatively large. For NxM pieces of puzzle we identify N points of fraction with the greatest value of D h function using [START_REF] Demaine | Jigsaw Puzzles, Edge Matching, and Polyomino Packing: Connections and Complexity[END_REF]. We select one of these points of fraction randomly using the biased roulette wheel selection mechanism so the point of fraction with the greatest dissimilarity has the highest probability to be selected. When the point of fraction has been determined we follow the horizontal direction till the next point of fraction is identified and the sequence of pieces between these two points of fraction constitutes the segment. An example of distribution of points of fraction as well as the relevant segments is on Fig. 2.

The final phase of this procedure is an attempt to allocate the segment into a new position within the picture using D h function value [START_REF] Demaine | Jigsaw Puzzles, Edge Matching, and Polyomino Packing: Connections and Complexity[END_REF] for the very first or very last piece of the segment. Moreover, we must not forget that albeit the point of fraction cannot occur at the beginning or at the end of rows (it results from the definition of D h ), the right place of segment could be at the beginning of the row and such a placement will be indicated by D v according to (2)). Of course, there is no guarantee that the sequence selected is a real and correct segment of the picture. However, as the points of fraction present very likely some faults within current arrangement of pieces, there is a chance to improve the overall composition.

It is clear that the above described asexual operator is efficient especially in latter phases of evolution when suitable building blocks (segments) exist within relevant individuals. That is why we have started the evolution process with 1000 individuals and carried it on for 50 generations. Then we employed decimation operation and discarded most of the population. Decimation is a secondary genetic operation that is normally performed at the beginning of a run and it is used to increase the proportion of fitter individuals in the population. It has been shown [START_REF] Nanduri | Comparison of the effectiveness of decimation and automatically defined functions[END_REF] that decimation provides improvement in performance that justifies the additional computation at the beginning of the run. Hence, taking advantage of the decimation operation only 100 individuals survived (we made sure that these individuals are unique in order to provide enough diversity) and then we start using the newly devised operator.

Experiments

We have experimented with three different pictures that were divided into NxN pieces. The size of the pieces was 80x80 pixels and so the size of the picture changed accordingly the number of pieces taken into account. The algorithm run for 200 generations, the first 50 generations using PMX operator with p c =0.8 and simple mutation (changing mutually places of two randomly selected pieces) with p m =0.05. For the next 150 generations the fraction point operator was applied to 50% of population whereas frequency of PMX slowed down to p c =0.3. Taking into account that the best individuals (10%) are copied automatically to the new population, the size of it remains constant for the rest of the evolution process. We ran our algorithm 50 times for each image and the results are given in Table 1. We can see that correct solution has been found for all the cases that were tested and the efficiency of the algorithm clearly deteriorates with increasing number of pieces. We have not included puzzles with more than 625 pieces (25x25) into our results table as we were unable to get a correct result within the given size of population and limited number of generations.

Discussion and Limitations

Our results show that the method using points of fraction is promising. There are still some opportunities to improve the way of points of fraction location as well as the way of matching the individual segments. Nevertheless, using this approach we were able to assemble puzzles of larger sizes than our predecessors.

In several cases our algorithm got stuck when it reached a partial solution in a form of horizontally cut picture that is glued together in a completely wrong way. We have analyzed this malfunction and we have realized that the algorithm was simply misled by a low level of dissimilarity along the cutting line. Because our algorithm does not make any use of the correct order and/or location of the individual pieces within the original picture, there is no way to avoid this kind of fault.

Conclusions

We have described the new algorithm for solving the rectangle piece color puzzle problem. Genetic algorithm has been augmented by new heuristic asexual operator that aims at identification of points of fraction within partially assembled picture, extraction of supposed sequence of correctly joint pieces, and its insertion into a new position in such a way that, if possible, the segment is enlarged. This method promotes the gradual enlargement of correctly connected blocks and the search for an optimal position of the block within the frame of the assembled picture is enabled and supported at the same time too. Our approach has been successfully tested and the algorithm is capable of solving puzzles consisting of several hundred pieces.

Fig. 1 .

 1 Fig. 1. Example of the puzzle used for the algorithm testing (10x10-piece puzzle).

Fig. 2 .

 2 Fig. 2. Points of fraction and prospective segments (10x10-piece puzzle).

Table 1 .

 1 Experimental Results.

	No. of pieces	Test 1	Test 2	Test 3
	15x15	100%	94%	100%
	20x20	92%	74%	86%
	25x25	56%	42%	52%


Ilias Bougoudis
email: ibougoudis@yahoo.gr

Lazaros Iliadis
email: liliadis@fmenr.duth.gr

Stephanos Spartalis

Comparison of Self Organizing Maps Clustering with Supervised Classification for air pollution data sets

Keywords: Self Organizing Maps, Artificial Neural Networks, Classification, Air Pollution

Air pollution is a serious problem of modern urban centers. The objective of this research is to investigate the problem by using Machine Learning techniques. It comprises of two parts. Firstly, it applies a well established Unsupervised Machine Learning approach (UML) namely Self Organizing Maps (SOM) for the clustering of Attica air quality big data vectors. This is done by using the concentrations of air pollutants (specific for each area) for a period of 13-years (2000-2012). Secondly, it employs a Supervised Machine Learning methodology (SML) by using multi layer Artificial Neural Networks (ML-ANN) to classify the same cases. Actually, the ANN models are used to evaluate the SOM reliability. This is done, because there is no actual and well accepted clustering of the related data to compare with the outcome of the SOM and this adds innovation merit to this paper.

Introduction

Air pollution has been defined as the release of any substances affecting the normal cycle of any vital process or degrade infrastructure [START_REF] Paschalidou | Neural Modeling of the Tropospheric Ozone concentrations in an Urban Site[END_REF]. The target of this research effort is to reveal existing air pollution patterns by performing distinct clusterings for each selected measuring station of the Attica area. The clusterings are based on a vast volume of air quality data vectors related to specific air pollutants for each site. The motivation to use Self Organizing Maps was mainly related to their potential to produce a low-dimensional (typically two-dimensional) discretized representation of the input space of the training samples. Another reason was the unsupervised learning mode of SOM [START_REF] Karatzas | Studying and predicting quality of life atmospheric parameters with the aid of computational intelligence methods[END_REF], [13]. It is well known, that the effectiveness of such an effort is obtained by estimating Sensitivity and Specificity indices as a result of the comparison between the actual clusters to the obtained ones [START_REF] Fawcelt | An Introduction to ROC Analysis[END_REF]. However, in the case examined here, there is no well accepted clustering either from the literature or from a study of the civil protection authorities. For this reason, Multi Layer Feed Forward (MLFF) ANNs were developed for each measuring location which classify the available data of each site. The input vectors of the MLFF approach comprised of air pollutants, plus seven meteorological and five daily factors. The outputs of the MLFF ANNs were considered as a comparison metric for the validation of the SOM's efficiency. It is the first time that SML is employed to support the validation of the UML classification. Herein, we will designate by classification the process of supervised learning on labeled data and by clustering the process of unsupervised learning on unlabeled data.

Literature review

There are several research efforts in the literature that use ANN to model pollution of the atmosphere [START_REF] Alkasassbeh | Prediction of PM10 and TSP Air Pollution Parameters Using Artificial Neural Network Autoregressive, External Input Models: A Case Study in Salt, Jordan[END_REF], [START_REF] Jordan | Neural Networks[END_REF]. [START_REF] Skön | Analysing Events and Anomalies in Indoor Air Quality Using Self-Organizing Maps[END_REF]. However, only recently some researchers have used SOM clustering to analyze air pollution of major cities. Patterns of air quality have been searched for Mexico City by Neme and Hernandez [START_REF] Paschalidou | Neural Modeling of the Tropospheric Ozone concentrations in an Urban Site[END_REF], whereas Karatzas and Voukantsis have done the same for the city of Thessaloniki [START_REF] Kirt | A method for comparing self-organizing maps: case studies of banking and linguistic data[END_REF]. Skön et al., have analyzed indoor air quality using SOM [START_REF] Yin | Learning Nonlinear Principal Manifolds by Self-Organising Maps[END_REF]. Li and Chou have investigated air pollution spatial variation with SOM [9]. Glorenec applied SOM to forecast Ozone peaks [START_REF] Glorennec | Forecasting Ozone Peaks Using Self-organizing Maps and Fuzzy Logic[END_REF]. Unfortunately though, due to different pollution measurements and different approach methods, any comparison of these works is inaccurate. Moreover, due to station's malfunctions in our case, there were a lot of records missing, which made our effort even more challenging. Also, we had only daily measurements for PM pollutants, while daily for every other pollutant. To the best of our knowledge there is no similar research for the wider area of Attica.

Materials and Methods

Area of research and Data

The area of research included four characteristic air pollution and meteorological stations. They were chosen as follows: "Athinas" is located exactly in the heart of the city centre. "Piraeus" has a unique position by the seaside, "Peristeri" is an urban site at Northwest and "Agia Paraskevi" is a suburb in the Northeastern part of Attica, close to a more mountainous area. So each measuring station is characterized by different topographic and geographic attributes. Figure 1 depicts the four locations under study. The available data sets that were used were related to the wider area of Athens, for a period of 13-years (2000-2012) and they varied from station to station. Overall, the pollutants of interest were hourly concentrations of monoxides plus dioxides like CO, NO , NO 2 , SO Wm ), wind speed ( sec m ) and direction (rad), pressure (mbar), Illumine and relative humidity. Finally, for every record we added its daily attributes; Year, Month, Day, Day_ID (1 for Monday, 7 for Sunday) and Hour. However not all stations measure the same pollutants. Thus, the inputs fed to both SOMs and MLFF ANNs were dependent on the sensors of each station. The meteorological parameters were the same everywhere and obtained from "Penteli" station, except from "Agia Paraskevi" station, where, because it is located between "Penteli" and "Thiseion" stations, we averaged the meteorological values for it.

Fig. 1. The four measuring stations

Data were obtained from the Greek ministry of Environment [START_REF] Iliadis | Artificial Neural Network Modelling of the surface Ozone concentration[END_REF]. Totally 1,017,733 vectors without missing values were available, whereas an average as high as 18.82% of the data are missing with the station of Piraeus having the worst percentage equal to 33.67%.

The following table 1 presents the types of parameters measured in each station and percentage of missing values [14].

Unsupervised Learning

In UML we let the algorithm decide how to group samples into classes that share common properties. Some examples of unsupervised learning are Kohonen's Self Organizing Maps, K-Means Clustering and Neural Gas Networks.

Self Organizing Maps (SOMs) are a well established unsupervised ML approach, based on competitive learning. Their main advantage is their ability to isolate clusters in high dimensional spaces [9]. The Self Organizing Map as all ANN consists of neurons (nodes). A weight vector is assigned to each neuron. This vector is of the same dimension as the input data vectors. The amount of nodes is the number of the clusters that will be used to group the input data. The obtained Map is a NxN space, where the data are scattered and arranged. The number of neurons is set as the square of the map. Their function can be summarized in four steps [START_REF] Li | Multi-Resolution Spatio-temporal Data Mining for the Study of Air Pollutant Regionalization Proceedings of the 33 rd Hawaii Conf. on System Sciences[END_REF]:

Initialization: All of the connection weights of each cluster are initialized Competition: In this stage for each input pattern, the neurons compete to each other in order to "win" this input. The neuron which adapts its value closest to the input "wins". We can define the discriminant function to be the squared Euclidean distance between the input vector x and the weight vector w j for each neuron j as:

  2 1 ) (     D i ji i j W X X d ( 1
)
Cooperation: Here follows the creation of a neighbourhood located close to the previously winning neuron. In this way, the winning neuron creates a neighbourhood with other neurons, in order to cooperate with each other and win future inputs. If S ij is the lateral distance between neurons i and j on the grid of neurons, we define a topological neighbourhood T j,I(X) where I(x) is the index of the winning neuron.

] 2 exp[ 2 2 ) ( , ) ( ,  X I j X I j S T   (2)
Adaption: In this last stage, each neuron creates a neighbourhood or becomes a member of a neighbourhood and self -organizes, so that the feature map between inputs is formed. In practice, the appropriate weight update equation is:

) ( ), ( ). ( ) ( , ji i X I j ji W X t T t n W    (3) [8].
In every step, all neurons adapt their weights to the current input, but not as much as the winner neuron and its neighbourhood [START_REF] Li | Multi-Resolution Spatio-temporal Data Mining for the Study of Air Pollutant Regionalization Proceedings of the 33 rd Hawaii Conf. on System Sciences[END_REF]. In this way, each neighbourhood is suitable for certain values, and so the map is ordered and shaped.

SOM Clustering

In this study, the number of neurons was the only subject of experimentation, in the development of the feature maps. The main target was to isolate the extreme pollutant values of every station and then to record the meteorological and spatiotemporal characteristics for these particular cases. Through this research we have experimented with 2, 3 and 4 neurons. As a result, the obtained maps have 4, 9 and 16 clusters accordingly. Often, the selection of a small number of clusters may lead to loss of the pattern (if there is one). However in this effort, the numbers mentioned above have proven to be suitable. Larger numbers produced larger maps, separating our input data in a chaotic level. For every map produced, we kept the hits of each cluster, the cluster of each record, the neighbour weight distances and the weights of each input.

The

Pattern Recognition with MLFF ANNs

Classification is a form of Supervised Learning. After we obtained the Clusters for each station, we used the pattern recognition tool of MATLAB to classify the available data records by developing MLFF ANNs (three for each site plus a daily one for the stations were we had PM pollutants measures). The purpose was to figure out if obtained grouping was suitable. For the development of the optimal ANNs the data were divided in training, evaluation and testing sets by 70%-30% and 30% respectively. The evaluation metrics were mean square error and confusion matrices. The training function was trainscg (Scaled conjugate gradient backpropagation) and the number of hidden neurons 10, for all networks.

As it has already been mentioned, the input data for every station comprised of the measured pollutants plus seven meteorological parameters (as shown in table 1) and five daily parameters (year, month, day, day identification and hour). By this approach, evaluation became stricter, as we had 17 input parameters totally. As it is shown in the following Figures, the performance of the ANN is very accurate for the "Peristeri" station for the period 2000-2004. If a pollutant was missing values (for example in Agia Paraskeui, the SO 2 pollutant has value only from 2000 to 2005), we created a separate SOM for it. Also, if 9 neurons could not group the extreme values of a pollutant, we developed a second SOM with more neurons (16), in order to manage to cluster the specific pollutant. Unfortunately, the confusion matrix of the networks with 16 clusters is not visible, as it contains 256 elements. For every network we obtained the confusion matrix and the percentages of correct and incorrect classifications. Figure 7 shows the confusion matrix which specifies the high compatibility level between clustering and classification for the "Peristeri" site. Unfortunately there is not enough space for all confusion matrices. However the high agreement of the two approaches is shown in the following tables.

Comparative analysis between SOMs and MLFF ANNs

The tables below show the results from both the clustering and the classification procedures. Columns 2 to 4 are the average, minimum and maximum values for each pollutant. Columns 5 to 8 are the same values extracted from each cluster which contained the extreme values for each pollutant. The last column is the correct percentage of the classification process.

Conclusions -Discussion

In almost every station we managed to isolate the extreme pollutants' values. From the clusterings we have obtained certain pollution patterns comprising of specific temporal, meteorological and air pollutant concentrations. While observing these conditions, we came to the conclusion that high values of CO and NO are present, when we have low temperature, high humidity, low solar radiation and low wind speed, between 8-12 AM or 6-10 PM. On the other hand, high levels of Ο 3 appear under high temperature, low humidity and high solar radiation and between 12 -6 PM.

In almost every case, the cluster which included the extreme values was in one corner of the map (in most cases, in the lower right cluster). This was because this cluster was isolated and not part of a neighbourhood.

Although most cases had their extreme values put together in the same cluster, the values of Ο 3 were separated from the others. What is more, the cluster which contained the extreme values of Ο 3 was across the one which contained the other pollutants.

Almost all percentages of correct classification were above 80%, showing that the SOM clustering can be considered reliable and compatible to the MLFF ANN classification.

In most cases, when we had two SOMS (with 9 and 16 clusters) the correct percentage of the clustering with 9 clusters is higher from the one with 16 clusters.

Most of the times, the network failed to assign the extreme records in the right neuron. Instead, it put every record in the closest cluster. Although this may sound like a malfunction of the classification, it may be used as a specific pattern recognition tool. This occurs because meteorological parameters act in a catalytic manner and they are not the primary source of pollutants.

This research is quite innovative and it has proven the ability of SOM to reveal clusters with specific attributes in real world problems and moreover their potential contribution for quality of life.

Future work will include the comparison and evaluation of SOM with other unsupervised methods (like fuzzy k-means and Neural Gas Networks), for the same case.

In the matrices below we have the daily attributes modes for every station's extreme pollutant values. Inside the parenthesis we have the counter for each mode (for example, most extreme records for NO 2 appeared in 2009 511 times, while they occurred on Friday (DAY_ID is 5)). PM 10 and PM 2.5 have given the same results, so we combined them into one column.

  following images show the results of the SOM clustering. They can represent quite efficiently almost every other case. The Neighbour Weight Distance Figure provides information about the neighbourhoods created; the darker colours represent larger distances, and the lighter colours represent smaller distances. In fact we have separated the neurons which were isolated and not part of any neighbourhood. These neurons had the extreme pollutant values and thus they correspond to very interesting cases.The Input Weights Figure shows the weight assigned to each node from every input. It must be specified that input1 to input5 correspond to CO, ΝΟ, ΝΟ 2 , Ο 3 , SO 2 respectively. Lighter colour means larger value. Here, wherever we have bright colours (yellow or orange), we have the isolated neurons who are assigned the extreme values. Another clue for the most extreme neuron was that this neuron usually had less record than the others (Figure4). The combination of these three figures helped us to figure out where the most hazardous values were. We are seeking yellow or orange from the Input's Weights Figure, while dark connections in the Distances Figure.The neuron in the bottom -left corner is the first, the one next to it is the second, while the last one is on the top -right corner.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Neighbor Distances for "Athinas" Hourly pollutants values, period 00-04 (the 1 st neuron seems to be the extreme one)

Fig. 5 .Fig. 6 .

 56 Fig. 5. Error Histogram for the MLFF ANN for the "Peristeri" station (00-04)

Fig. 7 .

 7 Fig. 7. Overall Confusion Matrix for the MLFF ANN for the "Peristeri" station (00-04)

Table 1 .

 1 Description of the stations employed for

this research ID Station's name Code Missing values Correct Data Vectors Station's data

	1	Ag. Paraskevi	AGP	12.32%	99,936	Ο 3 , ΝΟ, ΝΟ 2 , SΟ 2
	2	Athinas	ATH	21.86%	89,058	Ο 3 , ΝΟ, ΝΟ 2 , CO, SO 2
	3	Peristeri	PER	33.61%	75,668	Ο 3 , ΝΟ, ΝΟ 2 , CO, SO 2
	4	Piraeus	PIR	33.67%	75,600	Ο 3 ,ΝΟ, ΝΟ 2 , CO, SO 2
	5	Penteli	PEN	3.66%	109806	Meteorological
	6	Thiseion	THI	0.30%	113,632	Meteorological

Table 2 .

 2 Comparative analysis for "Agia Paraskevi" Station

	AGP	ALL RECORDS		SOM		CLASS
		AVG MIN MAX	AVG	MIN	MAX RMSE	%
	NO (00-04)	2,75	1	411	25,84	1	411
	NO 2 (00-04) O 3 (00-04)	20,15 91,38	1 1	347 330	95,25 190,67	45 155	347 330	0,1	93,6
	SO 2 (00-04)	7,35	2	244	25,98	2	244
	NO (05-08) NO 2 (05-08)	2,19 22,37	1 2	141 198	47,96 73,44	22 36	141 165	0,14	82,9
	O 3 (05-08)	76,34	1	279	166,9	143	279
	SO 2 (05-08)	6,02	2	100	32,2	24	100	0,06	98,1
	NO (09-12)	2,46	1	111	11,71	1	111
	NO 2 (09-12)	13,23	1	287	49,94	22	287	0,1	92,5
	O 3 (09-12)	86,91	2	251	146,01	129	251
	PM 10 PM 2,5	26,13 17,17	6 4	292 74	127,75 45,87	91 23	292	0,19	87,9

Table 3 .

 3 Comparative analysis for "Athinas" Station

	ATH	ALL RECORDS		SOM		CLASS
		AVG MIN MAX	AVG	MIN	MAX RMSE	%
	CO (00-04)	2,43	0,1	21,4	11,03	5	21,4	0,1	95,0
	NO (00-04)	72,86	1	908	659,62	545	908
	NO 2 (00-04)	73,83	1	377	184,14	128	377
	O 3 (00-04)	33,22	1	253	124,51	93	253	0,1	82,5
	SO 2 (00-04)	13,29	2	259	41,13	3	259
	CO (05-08)	1,78	0,2	11,8	7,93	4,1	11,8
	NO (05-08) O 3 (05-08)	57,57 30,97	1 1	787 199	536 -	451 -	787 -	0,09	95,8
	SO 2 (05-08)	8,76	2	126	34,55	2	122
	NO 2 (05-08)	63,92	4	275	147,13	117	275	0,15	43,8
	CO (09-12)	1,39	0,1	10,4	5,81	1,7	10,4
	NO (09-12)	50,72	1	678	431,76	362	678
	NO 2 (09-12)	55,35	3	323	94,4	39	263	0,09	97,5
	O 3 (09-12)	34,33	1	186	-	-	-
	SO 2 (09-12)	7,63	2	86	19,85	8	60

Table 4 .

 4 Comparative analysis for "Piraeus" Station

	PIR	ALL RECORDS		SOM	CLASS
		AVG MIN MAX	AVG	MIN	MAX RMSE	%
	CO (00-04)	1,58

Table 5 .

 5 Comparative analysis for "Peristeri" Station

	,12	90,3

Table 6 .

 6 Daily attributes for extreme pollutant values for "Agia Paraskevi" Station

	AGP	CO	NO	NO 2	O 3	SO 2	PM
	YEAR	2009 (511)	2009 (511)	2009 (511)	2010 (453)	2003 (247)	2008 (5)
	MONTH	1(253)	1(253)	1 (253)	8(499)	5(147)	4(5)
	DAY	16(74)	16(74)	16(74)	6(81)	30(47)	21(2)
	DAY_ID	5(210)	5(210)	5(210)	6(277)	4(160)	3(4)
	HOUR	8(159)	8(159)	8(159)	14(291)	10(127)	-

Table 7 .

 7 Daily attributes for extreme pollutant values for "Athinas" Station

	ATH	CO	NO	NO 2	O 3	SO 2	PM
	YEAR	2009 (139)	2009 (139)	2009 (139)	2002 (444)	2009 (139)	-
	MONTH	1(194)	1(194)	1(101)	7(363)	1(194)	-
	DAY	3(28)	3(28)	3(28)	4(61)	3(28)	-
	DAY_ID	5(72)	5(72)	5(72)	7(184)	5(72)	-
	HOUR	9(298)	9(298)	11(165)	16(154)	11(90)	-

Table 8 .

 8 Daily attributes for extreme pollutant values for "Piraeus" Station

	PIR	CO	NO	NO 2	O 3	SO 2	PM
	YEAR	2002 (145)	2002 (145)	2002 (294)	2002 (893)	2002 (294)	2007(11)
	MONTH	12(76)	12(76)	3(151)	5(294)	3(151)	12(7)
	DAY	30(27)	30(27)	26(36)	17(73)	26(36)	1(4)
	DAY_ID	3(59)	3(59)	1(95)	6(500)	1(95)	2(9)
	HOUR	9(85)	9(85)	10(136)	17(577)	10(136)	-

Table 9 .

 9 Daily attributes for extreme pollutant values for "Peristeri" Station

	PER	CO	NO	NO 2	O 3	SO 2	PM
	YEAR	2011 (282)	2011 (282)	2011 (436)	2003 (479)	2011 (282)	-
	MONTH	12(203)	12(203)	5(535)	7(579)	12(330)	-
	DAY	2(43)	2(43)	13(71)	24(61)	20(44)	-
	DAY_ID	5(113)	5(113)	5(231)	5(196)	5(113)	-
	HOUR	9(210)	9(210)	8(177)	15(310)	11(103)	-


Jonne Pohjankukka

Paavo Nevalainen

Tapio Pahikkala
email: tapio.pahikkala@utu.fi

Eija Hyvönen

Pekka Hänninen

Raimo Sutinen
email: raimo.sutinen@gtk.fi

Jari Ala-Ilomäki
email: jari.ala-ilomaki@metla.fi

Jukka Heikkonen
email: jukka.heikkonen@utu.fi

Predicting water permeability of the soil based on open data

Keywords: load bearing capacity of soil, water permeability, regression, k-nearest neighbor, mobility, sub-arctic infrastructure

Water permeability is a key concept when estimating load bearing capacity, mobility and infrastructure potential of a terrain. Northern sub-arctic areas have rather similar dominant soil types and thus prediction methods successful at Northern Finland may generalize to other arctic areas. In this paper we have predicted water permeability using publicly available natural resource data with regression analysis. The data categories used for regression were: airborne electro-magnetic and radiation, topographic height, national forest inventory data, and peat bog thickness. Various additional features were derived from original data to enable better predictions. The regression performances indicate that the prediction capability exists up to 120 meters from the closest direct measurement points. The results were measured using leave-one-out cross-validation with a dead zone between the training and testing data sets.

Introduction

This paper is about predicting the water permeability of the soil by regression analysis using publicly available multi-source data. Water permeability (also called hydraulic conductivity) is a central soil property related to soil type and soil texture. High permeability means that soil tends to stay dry and traversable most of the year, whereas low permeability creates a risk for mobility when precipitation is high. Mobility in arctic areas is of great interest to many different parties. E.g. the mining industry is interested about the mobility estimates when placing various facilities. The forest industry is interested on the load bearing capacity of the soil, since the route solutions can be adaptive to mobility predictions.

Our input data set consists of 44 features which are publicly available. The data is in raster format with grid resolution ranging from 10 meters to 50 meters.

Water permeability of the soil has been measured in 1788 test spots at Northern Finland provided by the Geological Survey of Finland (GTK). It is an important attribute which, when combined with other features available, helps to determine the soil types. Related studies have been conducted in [START_REF] Azzalini | Prediction of soil respiration rates from temperature, moisture and soil type[END_REF] where soil respiration rates are predicted from temperature, moisture content and soil type. Another related research was published in the paper of P. Scull, J. Franklin and O.A. Chadwick [START_REF] Sculla | The application of classification tree analysis to soil type prediction in a desert landscape[END_REF]. In their paper they use classification tree analysis for predicting the soil type in desert landscapes. R.P.O. Schulte et al. [START_REF] Schulte | Predicting the soil moisture conditions of irish grasslands[END_REF] focuses on soil moisture deficit, which is a related concept but not of concern in sub-polar areas, H. Gao et al. [START_REF] Gao | Algorithm Theoretical Basis Document for Terrestrial Water Cycle Data Records[END_REF] and R.A. Chapuis [START_REF] Chapuis | Predicting the saturated hydraulic conductivity of soils: a review[END_REF] focus on water budget modeling, which was not yet attempted in our study. H.S. Mahmood et al. [START_REF] Mahmood | Proximal gamma-ray spectroscopy to predict soil properties[END_REF] uses on-site gamma-ray measurements for analysis of the farming soil. N.J. McKenzie [START_REF] Mckenzie | Spatial prediction of soil properties using environmental correlation[END_REF] combines gamma-ray and digital elevation model to predict the chemical composition of the farming land. Closest to our paper is [START_REF] Emadi | An approach for land suitability evaluation using geostatistics, remote sensing, and geographic information system in arid and semiarid ecosystems[END_REF], where several data sources (topographic and remote sensing) are combined with 85 soil samples to assess the usability of the soil within and outside the sampled area. One can try to by-pass the water permeability estimation by directly learning the dynamic coupling between the precipitation and remotely observed soil moisture. This approach must include the digital terrain model (DTM) to estimate the water catchment. An example of this approach is [START_REF]Soil moisture estimation using remote sensing[END_REF].

The main novelty of this paper related to the previous studies is that the prediction is based on wide-area public data on a subpolar region. The features used in this paper are basically available through-out the arctic zone.

We use regression analysis to find a mapping between the publicly available data and water permeability of the soil. In the following, we present the regression methods in Ch. 2. Then we introduce the test area, the original data sets and derived features (Ch. 3) and describe the analysis process and results of the analysis (Ch. 4). The last part is for conclusions and future approaches (Ch. 5).

Regression methods

Regularized least squares regression (RLS) is well known so we describe it mainly to introduce the variables and the notation used later in the paper. The explanatory variables x 1 , ..., x p consist of given data and dependent variable y is the water permeability. We need to find a set of parameters w ∈ R p and b ∈ R such that the error function:

E(w) = 1 n n i=1 y i -w T x i -b 2 + λ n w T w (1)
is minimized, where x i ∈ R p is the input vector, y i ∈ R is the response value, n is the number of observations and λ is the regularization parameter.

The k-nearest neighbors (k-NN) approach predicts the test sample by taking the average from k points nearest to it. Euclidean distance is the used metric in our analysis. Explicitly stated, if y 1 , ..., y k are the response values of the knearest points to the test sample, then the response value for the test sample ŷt is:

ŷt = 1 k k i=1 y i .
3 Test area, data sets and features

The research area is located in the northern part of the municipality of Sodankylä, which is a part of Finnish Lapland. The size of the target area is 18432 km 2 . The center point of the rectangular target area is at ETRS-TM35FIN coordinates 7524 kmN, 488 kmE, zone 35.

The data set consists of aerial gamma-ray spectroscopy data (referred later as gamma-ray data, AGR) combined with electromagnetic (AEM), topographical (Z), peat bog mask (PBM) and The National Forest Inventory 2011 (VMI1 ) data when predicting the qualities and characteristics of the soil, namely its type and water permeability (WP). Gamma-ray data is inversely related on the amount of water on the soil, which can be used to predict the type of the soil. The forest inventory data describes the profile of tree species, their maturity and foresting state. Albeit this kind of data is not directly available elsewhere in northern subarctic areas (e.g. Russia, Canada), several studies are underway to predict the main characteristics of the forest by remote measurement methods [START_REF] Leutner | Modelling biodiversity and forest structure using hyperspectral[END_REF]. These methods include LiDAR and various satellite measurements.

The data providers are: When considering all the derived features used in the analysis we get a total of 96 data layers.

The test site has 1788 sample points, where many mechanical and electrochemical properties of the soil were measured, see [START_REF] Hyvönen | Assessing site suitability for scots pine using airborne and terrestrial gamma-ray measurements in finnish lapland[END_REF]. The water permeability is a theoretical value derived from the soil particle size distribution of the soil.

We now present our data sources and donors.

Forest inventory data

The National Forest Inventory (VMI) holds the state of Finnish forests. The data is updated once in two years. The parameters are derived from various remote sensing sources, and several spot-wise verification and calibration methods are applied to it before publishing the data [START_REF] Tomppo | Multi-source national forest inventory methods and applications[END_REF]. 44 numerical features include green mass, trunk dimensions and tree density per specie category. These multi-source features exhibit built-in dependencies, thus the final number of useful features is lower.

Aerial gamma-ray data

The aerial gamma-ray data was provided by the Geological Survey of Finland (GTK). The raster data is based on gamma-ray flux from potassium, which is the decay process of the naturally occurring chemical element potassium (K). This data indicates many significant characteristics of the soil, including the tendency to stay moist after precipitation and tendency to frost heaving. Also the soil type, especially density, porosity, grain size and humidity of the soil have an effect to gamma-ray radiation. In Fig. 1 we present the gamma-ray data from Sodankylä target area. The bright end of the gray scale is for the high gamma radiation and hence less water in the locality of the pixel.

Electromagnetic properties of soil

The air-borne electromagnetic (AEM) data was provided by the Geological Survey of Finland (GTK). Primary AEM components, in-phase and quadrature, were transformed to apparent resistivity values by using a half-space model [START_REF] Hautaniemi | The three in one aerogeophysical concept of gtk in 2004[END_REF].

The apparent resistivity gives information on different kind of soil conductors. The apparent resistivity is governed by grain size distribution, water and electronic conductors content of soil and cumulative weathering.

Topographical height data

Topographical data provided by the National Land Survey of Finland (NLS) was included in the analysis. The data from NLS server is basically similar to aerial laser measurements (LiDAR) except LiDAR can reach denser grid. Instead of raw height alone we used local height difference, flow accumulation area, confluence and inclination described in [START_REF] Schwanghart | Topotoolbox: a set of matlab functions for topographic analysis[END_REF]. These four derived features are more efficient for prediction than raw height data alone.

Peat bog mask

Peat bog mask is created from GTK aero-radiometric data and is courtesy of NLS and METLA. The grid size is 20 m and the value 1 indicates that peat thickness is over 60 cm. Value 0 indicates thickness less than 60 cm. The limit chosen is practical for mobility prediction.

Derived features

The following features were derived from gamma-ray and electromagnetic data:

-Mean and variance over 3 × 3 window -Mean and variance over Gabor filter with 8 orientations, see [START_REF] Weldon | Efficient gabor filter design for texture segmentation[END_REF] -Local Binary Pattern (LBP) with pixel radii r ∈ {1, 2}, see [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture[END_REF] From topographical height we derived the following features: local height difference, ground inclination, convergence index and flow accumulation area. The definition of these features is at [START_REF] Kiss | Determination of drainage network in digital elevation models, utilities and limitations[END_REF]).

There are several additional attributes possible to derive from topographical height data, and more geomorphological features will be employed in the future.

The regression methods use total of 44 original and 52 derived features, including the constant feature. The derived features are useful only if the original feature is continuous enough. E.g. the Forest Inventory data often has locally constant zones with abrupt changes and the derived features do not help much.

Water permeability exponent

This is the subject of prediction. Basically, the water permeability indicates the nominal vertical speed of water through the soil sample. The measurement of this quantity is indirect, based on soil particle size distribution, and the actual speed highly depends on the inhomogeneities (roots, rocks) and micro-cracks in the soil. This is why this quantity is descriptive and theoretical. In our analysis we are using a logarithmic quantity x wp derived from water permeability speed v. For purposes of this presentation it is called as the water permeability exponent and defined as:

x wp = -log 10 v, [v] = m sec , (2)
This formula has v as the vertical speed of water flow through the soil.

Analysis and results

We are looking for methods which predict water permeability on areas, where there may not be direct water permeability measurements nearby. Therefore, we developed a modification of the leave-one-out cross-validation (LOOCV) for measuring the degree of spatial dependency from the nearby direct measurements, which we refer to as LOOCV with dead zone. Namely, the approach works on the measurement data just like an ordinary LOOCV in which each measurement at a time is omitted from the training set and used as a test point, except that we also remove from the training set all points that are within geographical distance r from the test point. This approach is illustrated in Fig. 2.

By varying r, we can measure how far from the test area we assume the closest measurements to be at the very least. In addition, the results can be helpful in deciding how dense grid of direct measurement one should use in order to obtain a certain level of prediction performance. We perform the regression of water permeability with the following three feature sets:

-location only -features + location -features only where location refers to the geographical coordinates (e.g. latitude and longitude) and features to the ones described in Section 3. Note that one can not rely on the location information if there are no nearby direct measurements at all, and therefore we measure the prediction performance separately with these.

The prediction performances with the different feature sets as a function of the radius r of the dead zone are depicted in the two leftmost graphs in Fig. 3 on p. 8. The generic version based on feature data only gives weaker results, since the sample point arrangement at Sodankylä (see sample sets A and B in Fig. 2) and perhaps the phenomenon itself induce spatial dependency. No good generic regression method for this data set has been found, instead the problem is about how much additional samples are needed per target area to make the prediction useful.

The common k-NN method has one essential parameter, the number of neighbors k. The spatial dependency can be probed by adding the dead zone radius r to avoid the optimistic effect of the nearest neighbors. Fig. 2 depicts the modified leave-one-out arrangement, where k nearest points outside the dead zone of radius r are used for teaching. By varying r one gets a varied data set and a rough estimate on how dense it should be for it to predict well in new circumstances.

The same parameterized dead zone leave-one-out arrangement was used with regression, too.

Predicting water permeability

As mentioned in Sec. 3.7 before, the prediction subject is the water permeability exponent x wp defined by Eq. 2. The values used for regularization parameter λ ranged from 2 -15 , ..., 2 15 . k-NN parameter had k ∈ {1, 3, 6, 12, 22}. Two different error measures were used for estimating prediction performance: mean absolute error (MAE) and concordance index (CI) [START_REF] Gönen | Concordance probability and discriminatory power in proportional[END_REF]. Explicitly, the error measures are:

M AE = 1 n n i=1 y i -ŷi y i CI = 1 N y i <y j h(ŷ i -ŷj ) (3)
MAE prediction baseline ỹ is the best possible prediction under the assumption that the prediction will be constant, thus constraining all values ŷj = ỹ, j = 1..n in the error minimization process. MAE baseline becomes thus:

M AE b = arg min ỹ 1 n n i=1 | y i - ỹ y i |
The prediction performance should be better than this to be useful. The corresponding percentage values (MAPE and MAPE b ) have been used in the rest of the text.

Concordance index counts the occurrences when the prediction fails to be monotonical. In equation (3) we denote N = | {(i, j) | y i > y j } | as the normalization constant which equals to the number of data pairs with different label values. h(.) is Heaviside step function.

The values of the C-index range between 0.0 and 1.0, where 0.5 corresponds to a random predictor and 1.0 to the case where prediction is monotonically correct.

Results

The results for regression analysis can be seen in Fig. 3 on p. 8.

Both MAPE and C-index indicate rather good prediction performance to the distance of 120 m from the nearest soil sample point. This is seen both with k-NN and RLS methods. When MAPE is higher than the baseline, it is better to use baseline average than the prediction. MAPE baseline is the horizontal line in the lower figures in Fig. 3. The dead zone radius r > 0 simulates a situation, where the test point is at least r distance away from the given training points. r = 0 is traditional LOO test arrangement and measures best the properties of the predicted value within the training set itself. It may be too optimistic, since we seek for generalization. A large radius r ≈ ∞ is overly pessimistic, since it would use only tiny fragments of the training set and would completely distort the prediction.

The prediction performance near r = 0 seems to indicate rather good generalization ability, but the performance reduces drastically over the dead zone distance r. Further study, both theoretical and practical, must be done to prop-erly address classical geoinformatics concepts such as spatial autocorrelation and spatial semivariance together with the general prediction ability. The problem is new, since spatial analysis in geosciences is usually applied in the sense of interpolation and extrapolation performance, and general prediction is usually analyzed in the terms of Machine Learning performance.

The feature selection was not attempted. There were two reasons for this:

-the number of features (96) remained modest.

-the forest inventory features are unique to Finland. They can be largely substituted by various remote measurements [START_REF] Tomppo | Multi-source national forest inventory methods and applications[END_REF], which would extend the application scope of the method to whole sub-polar area.

Conclusions and future work

The results indicate that the chosen five data sources (forest inventory, gammaray, air-borne electromagnetic, topographical data and peat bog mask) can be used to estimate the water permeability to a certain range from known measurements. This range seems to be c. 120-150 m. The best results come from the k-NN method based on the location of the sample points only. This method is naturally unavailable for general prediction. There are several possible improvements. Since the mapping from water permeability to soil types is not unique, see [START_REF] Pohjankukka | Arctic soil hydraulic conductivity and soil type recognition based on aerial gamma-ray spectroscopy and topographical data[END_REF], a special majority rule could be used to select the dominant soil type from neighboring grid point predictions. Such expert rules would require additional features like sophisticated geomorphological categories.

The Aerial Light Detection and Ranging (LiDAR) data can substitute most of the topographical and forest inventory data features. This would extend the scope of the prediction to any location at the arctic zone, where only aerial and satellite measurements are economical. LiDAR has also potential for derived features like geological morphology [START_REF] Sutinen | Airborne lidar detection of postglacial faults and pulju moraine[END_REF] and soil water budget modeling [START_REF] Leutner | Modelling biodiversity and forest structure using hyperspectral[END_REF].

The final goal is to predict the water permeability, soil types, approximate water budget and the load bearing capacity of the terrain in relation to the given weather forecast, while the model is based on remote measures and online learning based on measurements from the harvester fleet. The potential applications aim to wide-area routing and location planning. In this regard, even a modest prediction power of features-only prediction could yield a cumulative effect on route decisions.

Fig. 1 :

 1 Fig.1: Aerial data: gamma-ray (left) and electromagnetic data (right). Air-borne electromagnetic data is sensitive to geological properties to depth of hundreds of meters, but it also indicates some features of the top soil.

Fig. 2 :

 2 Fig. 2: Left: 1788 sample points. Set A (1187 points, marked with red circles, distance to the nearest neighbor dNN ≤ 86m) is tightly packed and set B is very sparse (601 points, marked with black dots, aver. dNN ≈ 1.1 km). Right : the dead zone (with radius r) around the leave-one-out test point (black circle). The gray circles are omitted from the training set (white circles). Both the k-NN and RLS method address the training data only, e.g. the k nearest neighbours are selected from outside the circle.

Fig. 3 :

 3 Fig. 3: Left : k-NN results with k = 6 and 3 different feature sets. Right : RLS results on features-only case. C-index and Pearson correlation at top and MAPE below. The prediction performance is adequate below 120 meters.

Table 1 :

 1 Data providers, data and the grid size.

	Provider	Data	Grid size
	Geological Survey of Finland (GTK) AGR, AEM	50 m
		WP
	Finnish Forest Research Institute (Metla) VMI, PBM	20 m
	National Land Survey of Finland (NLS)	Z	10 m

VMI2011: http://www.metla.fi/ohjelma/vmi/vm11-info-en.html

ACKNOWLEDGEMENTS

This work is done as a part of ULJATH project, which is funded by the Finnish Funding Agency for Technology and Innovation (TEKES). ULJATH stands for New Computational Methods For The Efficient Utilization of Public Data.

The authors would like to thank the anonymous reviewers for valuable comments and suggestions.