Van Do

Keisuke Nguyen
email: va-nguyen@kddilabs.jp

Kazunori Ogawa

Masayuki Matsumoto

Hashimoto

Editing Training Sets from Imbalanced Data Using Fuzzy-Rough Sets

Keywords: Rough Set theory, fuzzy-rough sets, classification performance, instance selection, imbalanced data

In this research, we study several instance selection methods based on rough set theory and propose an approach able to deal with inconsistency caused by noise and imbalanced data. Recent attention has focused on the significant results obtained in selecting instances from noisy data using fuzzy-rough sets. For imbalanced data, fuzzy-rough sets approach is also applied before and after using balancing methods in order to improve classification performance. In this study, we propose an approach that uses different criteria for minority and majority classes in fuzzy-rough instance selection. It thus eliminates the step of using balancing techniques employed in controversial approach. We also carry out some experiments, measure classification performance and make comparisons with other methods.

Introduction

Rough set model is a machine learning technique for acquiring knowledge from datasets. Rough set theory was first introduced by Pawlak [START_REF] Pawlak | Rough sets[END_REF][START_REF] Pawlak | Rough Sets. Theoretical Aspects of Reasoning about Data[END_REF] in the early 1980s. It is a mathematical approach designed to deal with the ambiguity, vagueness and uncertainty that exist in a database. In knowledge discovery, rough set theory constitutes a sound basis by offering mathematical tools to discover patterns hidden in data. It can, according to [START_REF] Skowron | The discernibility matrices and functions in information systems[END_REF], be used for feature selection, data reduction, decision rule generation, and pattern extraction.

Recently, rough sets have been implemented for instance selection [START_REF] Caballero | Improving the k-nn method: Rough set in edit training set[END_REF][START_REF] Jensen | Fuzzy-rough instance selection[END_REF][START_REF] Verbiest | Frps: A fuzzy rough prototype selection method[END_REF]. The main ideas of these approaches are extracting fuzzy memberships of positive regions and choosing the instances which have large membership degrees for training phases. However, it is believed that these criteria do not overcome all problems in relation to inconsistency issues.

One issue that may cause an inconsistency problem is imbalanced data [START_REF] Lopez | An insight into classification with imbalanced data: Empirical results and current trends on using data intrinsic characteristics[END_REF]. A dataset is imbalanced when the numbers of instances in some classes are much larger than in others. Such classes are called majority classes. The classes with small cardinality are referred to as minority classes.

In fact, some researchers are trying to use the advantages of rough sets to enhance classification performance when a balancing method is employed [START_REF] Ramentol | *: a hybrid preprocessing approach based on oversampling and undersampling for high imbalanced data-sets using SMOTE and rough sets theory[END_REF][START_REF] Ramentol | Smote-frst: a new resampling method using fuzzy rough set theory[END_REF][START_REF] Verbiest | Preprocessing noisy imbalanced datasets using SMOTE enhanced with fuzzy rough prototype selection[END_REF]. In these approaches, fuzzy-rough sets are first used to remove low quality instances. Then, a well-known balancing technique called "Synthetic Minority Oversampling Technique (SMOTE)" [START_REF] Chawla | Synthetic minority over-sampling technique[END_REF] is employed to form a candidate set. Finally fuzzy-rough sets are used again to select quality instances from the candidate sets. Obviously, such methods need considerable computation time to edit and select instances. This research will introduce an approach that uses fuzzy-rough sets to balance and select quality instances from imbalanced datasets. In this work, different thresholds for majority and minority classes are used, thus, more items from minority classes can be selected. The experiment results show considerable improvement in classification performance and the advantage with respect to other methods in the literature.

The paper is structured as follows: Section 2 reviews the original rough set theory and an extension called fuzzy-rough set. In section 3, some rough set approaches for dealing with inconsistency are discussed along with their issues. The new algorithms and experiments are introduced in Section 4. Finally, we present our conclusions in Section 5.

Theoretical Background

Information Systems and Equivalence Relation

An information system is represented as a data table. Each row of this table represents an instance of an object such as people, things, etc. Information about every object is described by object attribute (feature) values.

An information system in the rough sets study is formally defined as a pair I = (U, A), where U is a non-empty finite set of objects called the universe and A is a non-empty finite set of attributes such that f a : U → V a for every a ∈ A [START_REF] Pawlak | Rough sets[END_REF][START_REF] Pawlak | Rough Sets. Theoretical Aspects of Reasoning about Data[END_REF]. The non-empty discrete value set V a is called the domain of a. The original rough set theory deals with complete information systems in which ∀x ∈ U , a ∈ A, f a (x) is a precise value.

Any information system taking the form I = (U, A ∪ {d}) is called a decision table where d / ∈ A is called a decision (or label) and elements of A are called conditions. Let V d = {d 1 , ..., d k } denote the value set of the decision attribute, decision d then determines a set of partitions {C 1 , C 2 , ..., C k } of universe U , where

C i = {x ∈ U |f d (x) = d i }, 1 ≤ i ≤ k. Set C i is called the i -th decision class or concept on U . We assume that every object in U has a certain decision value in V d .
Formally, in complete information systems, relation EQU P (x, y), P ⊆ A, denotes a binary relation between objects that are equivalent in terms of values of attributes in P [START_REF] Pawlak | Rough sets[END_REF]. The equivalence relation is reflexive, symmetric, and transitive. Let E P (x) = {y ∈ U |EQU P (y, x)} be the set of all objects that are equivalent to x by P , which is then called an equivalence class. The family of all equivalence classes (or partitions) on U based on an equivalence relation is referred to as a category and is denoted by U/EQU P .

From equivalence classes, Pawlak [START_REF] Pawlak | Rough sets[END_REF][START_REF] Pawlak | Rough Sets. Theoretical Aspects of Reasoning about Data[END_REF] defined an approximation space that contains lower and upper approximations denoted by apprX and apprX, respectively, of set X ⊆ U as follows:

appr P X = {E P (x)|x ∈ U, E P (x) ⊆ X} = {x ∈ U |E P (x) ⊆ X}, (1)
appr P X = {E P (x)|x ∈ U, E P (x) ∩ X = ∅} = {x ∈ U |E P (x) ∩ X = ∅}. (2)
In decision table I = (U, A ∪ {d}), the positive region P OS P (X) of set X in terms of attribute set P ⊆ A is defined as:

P OS P (X) = {appr P E {d} (x)|x ∈ X}. (3)
Others region definitions such as negative, boundary and rough set properties can be obtained from the original rough set theory research [START_REF] Pawlak | Rough sets[END_REF][START_REF] Pawlak | Rough Sets. Theoretical Aspects of Reasoning about Data[END_REF].

Apart from using equivalence relations to define a rough set in complete information systems, there are also numerous studies [START_REF] Grzymala-Busse | Generalized probabilistic approximations of incomplete data[END_REF][START_REF] Kryszkiewicz | Rough set approach to incomplete information systems[END_REF][START_REF] Nguyen | On probability of matching in probabilistiy based rough set definitions[END_REF][START_REF] Nguyen | Rough set approach with imperfect data based on dempster-shafer theory[END_REF][START_REF] Nguyen | Discretization problem for rough sets methods[END_REF] that deal with incomplete or imperfect information systems in which data are not described by precise and crisp values.

Fuzzy Rough Set

A classical (crisp) set is normally defined as a collection of elements x ∈ X that can be finite, countable or over countable. Each single element can either belong to or not belong to a set X ⊆ X. For a fuzzy set, a characteristic function allows various degrees of membership for elements of a given set. Then a fuzzy set X in X is a set of ordered pairs [START_REF] Zimmermann | Fuzzy Set Theory and its Applications[END_REF]:

X = {(x, µ X (x))|x ∈ X}, (4)
where µ X (x) ∈ [0, 1] is called the membership function or grade of membership (also the degree of compatibility or degree of truth) of x in X . In crisp sets, for two special sets ∅ and U the approximations are simply defined as µ appr R U (x) = 1 and µ appr R ∅ (x) = 0. Based on the two equivalent definitions, lower and upper approximations may be interpreted as follows: An element x belongs to the lower approximation appr R X if all elements equivalent to x belong to X. In other words, x belongs to the lower approximation of X if any element not in X is not equivalent to x, namely, µ R (x, y) = 0. Likewise, x belongs to the upper approximation of X if µ R (x, y) = 1. Now the notion of rough sets in crisp sets is extended to included fuzzy sets. Let µ X and µ R denote the membership functions of set X and of the set {(x, y) ∈ U × U |R(x, y)}, respectively. The fuzzy approximation space [START_REF] Radzikowska | A comparative study of fuzzy rough sets[END_REF] of fuzzy set X on X in terms of fuzzy relation R(x, y) can be defined as follows:

µ appr R X (x) = inf I{µ R (x, y), µ X (y)}, (5)
µ appr R X (x) = sup T {µ R (x, y), µ X (y)}. ( 6
)
where I and T are fuzzy implicators1 and triangular norms (t-norm)2 , respectively.

From the above equations, it is noticed that membership of lower and upper approximation for one instance strongly depend on only one other instance since they are defined by minimum and maximum. To soften this definition, Cornelis et. al. [START_REF] Cornelis | Ordered weighted average based fuzzy rough sets[END_REF] suggest a fuzzy-rough set definition using order weight averaging (OWA) aggregation [START_REF] Yager | On ordered weighted averaging aggregation operators in multicriteria decisionmaking[END_REF]. An OWA operator F W is a mapping F : R n → R using weight vector W =< w 1 , ..., w n > such that

F W (a 1 , • • • , a n ) = n i=1 w i b i , (7)
where b i is the i th largest of the values in {a 1 , • • • , a n }.
An OWA operator is bounded by the minimum and maximum of vector W . It thus can soften minimum and maximum operators. Suppose that there are two vectors

W min =< w min 1 , • • • , w min n > and W max =< w max 1 , • • • , w max n >, where w min 1 ≤ • • • ≤ w min n and w max 1 ≥ • • • ≥ w max n
, we can define OWA fuzzy rough sets as follows:

µ appr R X (x) = F W min (I{µ R (x, y), µ X (y)}), (8)
µ appr R X (x) = F W max (T {µ R (x, y), µ X (y)}). (9)
In fact, the definition of a fuzzy relation on an attribute set depends on individual systems. One example of defining relation methods to deal with incomplete information systems is shown in [START_REF] Nguyen | Extended tolerance relation to define a new rough set model in incomplete information systems[END_REF]. Further investigation of the combination between fuzzy and rough sets can be found in [START_REF] Cornelis | Ordered weighted average based fuzzy rough sets[END_REF][START_REF] Dubois | Rough fuzzy sets and fuzzy rough sets[END_REF][START_REF] Dubois | Putting rough sets and fuzzy sets together[END_REF][START_REF] Radzikowska | A comparative study of fuzzy rough sets[END_REF][START_REF] Yao | Combination of rough and fuzzy sets based on -level sets[END_REF].

Rough Selection and Problem Stated

Significant research on instance selection using rough sets model was introduced the work of Caballero et al. [START_REF] Caballero | Improving the k-nn method: Rough set in edit training set[END_REF]. For this purpose, the authors calculated approximations and positive regions for each class in training sets. Instances for training phases are then selected in two ways. In the first method, they try to delete all instances in the boundary. The second employs a nearest neighbourhood algorithm to relabel instances in the boundary region. The issues associated with these approaches were raised in Jensen and Cornelis' study [START_REF] Jensen | Fuzzy-rough instance selection[END_REF]. Another possible limitation is that they just deal with crisp sets of attributes in decision tables.

To deal with the problems identified in the study of Caballero et al., an approach called Fuzzy-Rough Instance Selection (FRIS) has been proposed [START_REF] Jensen | Fuzzy-rough instance selection[END_REF]. The notion of this approach is using memberships of positive regions to determine which instances should be kept and which instances should be discarded for training. First, the method calculates relations among instances and then measures quality of each instance by its positive region membership. An instance will be removed if the membership degree is less than a specified threshold. For the same purpose, another approach of measuring quality of instances was introduced in [START_REF] Verbiest | Frps: A fuzzy rough prototype selection method[END_REF]. Such approaches will undoubtedly refine the positive regions with quality instances.

However, one limitation to be considered is that there is only one threshold used for all classes. In some applications such as medical disease prediction, it is not uncommon to be interested in some specific groups, disease groups for example, rather than the others (healthy groups). When only one threshold is used, once noise is eliminated, valuable instances in interesting classes also disappear.

Some approaches use the advantage of fuzzy-rough to qualify instance sets formed by SMOTE [START_REF] Ramentol | *: a hybrid preprocessing approach based on oversampling and undersampling for high imbalanced data-sets using SMOTE and rough sets theory[END_REF][START_REF] Ramentol | Smote-frst: a new resampling method using fuzzy rough set theory[END_REF][START_REF] Verbiest | Preprocessing noisy imbalanced datasets using SMOTE enhanced with fuzzy rough prototype selection[END_REF]. Therefore, the issues above may be avoided. The steps of these methodologies can be summarized as follows:

1 Using fuzzy-rough sets to calculate the qualities of every instance and choose instances having high quality. 2 Using SMOTE to create artificial instances from the first step and add to datasets 3 Using fuzzy-rough sets to eliminate low quality instances from datasets after the second step.

There are some differences among these methods. In Ramentol et al.'s study named SMOTE-RSB * [START_REF] Ramentol | *: a hybrid preprocessing approach based on oversampling and undersampling for high imbalanced data-sets using SMOTE and rough sets theory[END_REF][START_REF] Ramentol | Smote-frst: a new resampling method using fuzzy rough set theory[END_REF], the first step does not appear. In step 3, the algorithm removes only low quality instances in the artificial instance sets. In Verbiest's research [START_REF] Verbiest | Preprocessing noisy imbalanced datasets using SMOTE enhanced with fuzzy rough prototype selection[END_REF], quality measurements are different for step 1 and 3. The first measures performance for imbalanced data while the last deals with balanced information. Thus, they call the algorithm FRIPS-SMOTE-FRBPS (hereinafter referred to as FRPSS in this paper). It is also claimed in [START_REF] Verbiest | Preprocessing noisy imbalanced datasets using SMOTE enhanced with fuzzy rough prototype selection[END_REF] that FRPSS produces better training sets compared with SMOTE-RSB * .

It appears that FRPSS is an excellent tool for editing training sets. However, applying three instance selections steps may require considerable computation time compared with other methods.

Multiple Thresholds Fuzzy Rough Instance Selection

As stated in [START_REF] Jensen | Fuzzy-rough instance selection[END_REF], the Fuzzy-Rough Instance Selection (FRIS) algorithm is quite efficient as it eliminates some instances that might affect the positive region membership of the remaining objects. It thus refines the dataset by choosing quality instances for training. However, we also noted some limitations of conventional methods. An example is shown in Figure 1. If we use fuzzy lower approximation as instance qualities and compare these numbers with the threshold = 0.7 (FRIS method), numerous positive instances (minority class) will be removed while almost all negative instances will be remained. The same problem occurs with FRPSS in which an average of approximations is used as a quality measurement. Therefore, in this section, using the advantages of fuzzy rough set instance selection approach, we propose an instance selection method which uses different thresholds for different classes of decision tables. By means of this strategy, SMOTE or other balancing techniques may not be necessary.

MFRIS algorithms

In decision table I = (U, A ∪ {d}), let R a (x, y) and µ Ra (x, y) denote a fuzzy relation and membership function, respectively, between objects x, y ∈ U on attribute a ∈ A. Then the membership function of the relation on an attribute set P ⊆ A is defined as:

µ R P (x, y) = T a∈P µ Ra (x, y),
where T is a t-norm.

In fact, there are many ways to calculate the membership function of a fuzzy relation between two instances. It is also possible to use distance based similarity such as the normalization of Euclidean distance. It depends on the characteristics of each system.

Data:

Xi, Yj, minority and majority classes; tX , tY , the selection thresholds for minority and majority classes. Result: Decision table (S, A ∪ {d}) Calculate membership degree of the relation among objects: µR(x, y); Calculate membership degree of the approximations for each class; Calculate the quality measurement of all instances for their classes; S ← ∅; The membership functions of approximation spaces for an object set X ⊆ U on attribute P ⊆ A can be defined by either equations ( 5) and ( 6) or ( 8) and [START_REF] Dubois | Putting rough sets and fuzzy sets together[END_REF]. In this study, as mentioned in Section 2, we must note that in the decision table

for x ∈ Xi do if µ R X i (x) ≥ tX then S ← S ∪ {x}; end end for x ∈ Yj do if µ R Y i (x) ≥ tY then S ← S ∪ {x}; end end
I = (U, A ∪ {d}) is a single label decision table such that f d (x) is a precise value.
In fuzzy sets, for a single label decision table, µ Ci (x) > 0 if f d (x) = d i and µ Ci (x) = 0, otherwise.

In selecting learning instances, we first define a function to measure the quality of an instance for each class X in terms of relation R:

µ R X (x) : αµ appr R X (x) + (1 -α)µ appr R X (x), (10)
where α ∈ [0, 1], X is the fuzzy set on X. In the above equation, note that x may or may not belong to X. This is where this approach differs from other approaches. When comparing µ R X (x) for x ∈ X, we can assume that α = 1 as in the study of Jensen et al. [START_REF] Jensen | Fuzzy-rough instance selection[END_REF] and α = 0.5 as in the approach of Verbiest et al. [START_REF] Verbiest | Preprocessing noisy imbalanced datasets using SMOTE enhanced with fuzzy rough prototype selection[END_REF]. Now, let t X and t Y be the selection thresholds for minority and majority classes, X i and Y j be the fuzzy sets on minority class X i and majority classes Y j , respectively, R denotes the fuzzy relation on the universal. The set of instances selected for the training phase can be defined as:

S = ∪ i {x ∈ X i |µ R Xi (x) ≥ t X } ∪ ∪ j {x ∈ Y j |µ R Yi (x) ≥ t Y } . (11)
Then, the algorithm to select instances can be described as shown in Figure 2.

In the first algorithm, depending on labels of instances, the quality measurement of every instance on their classes will be compared with the threshold for minority or majority thresholds. This means that an instance will be deleted from a training set if it is low quality even it is in minority classes. The different thresholds may help us to keep more instance in minority classes if necessary.

In addition, we also propose the second algorithm to select instances. The notion of this method is keeping all instances in minority classes while removing or relabelling some instances in majority classes. To describe the algorithm, we first introduce some definitions.

Let X, Y be the family of minority and majority classes respectively. Let X i , Y j denote fuzzy sets on X i , Y j , respectively. Thresholds for minority and majority classes are t X and t Y , respectively. The set of instances for which labels could be changed can be defined as

S Y →X = i {x ∈ Y i |µ R Yi (x) < t Y ∧ sup j {µ R Xj (x) -t X } ≥ 0}. ( 12
)
From the above definition, there are some instances in majority classes whose labels could be changed to a minority class label. The possible minority class set can be determined as follows:

Z(x) = i {X i |µ R Xi (x) ≥ t X , x ∈ S Y →X }. (13)
Then we can re-calculate the class membership functions of an instance x ∈ S Y →X ∩ Y n as µ Xm (x) = µ Yn (x); µ Yn (x) = 0, where m is defined to satisfy:

X m ∈ Z(x); µ R Xm (x) ≥ µ R Xi (x), ∀i, X i ∈ Z(x)
. Finally, the selected instances for the training phase can be defined as:

S = ( i X i ) ∪ ( i {x ∈ Y i |µ R Yi (x) ≥ t Y }) ∪ S Y →X . ( 14
)
The algorithm is then described as shown in Figure 3.

Experiment

Experiments were conducted on datasets from KEEL dataset repository [START_REF] Alcala-Fdez | Keel data-mining software tool: Data set repository, integration of algorithms and experimental analysis framework[END_REF] and UCI Machine Learning Repository [START_REF] Bache | UCI Machine Learning Repository[END_REF] and our lifestyle-diseases dataset. The dataset lifestyle-diseases is raw data of collected health factors such as age, BMI, and waist circumference,... and the outcomes are related to lifestyle diseases. Datasets and their properties are shown in Table 1. In this table, "IR" shows imbalance rates for the number of instances between majority and minority classes.

Each dataset is divided into three parts for both minority and majority classes. For each train-test, we used two parts to form training data and the rest were used as testing instances.

Data: X, Y , the family set of minority and majority classes tX , tY , the selection thresholds for minority and majority classes Result: Decision table (S, A ∪ {d}) Calculate membership degree of the relation among objects: µR(x, y); Calculate membership degree of the approximations for each class; Calculate quality measurement of all instances in the majority classes for all classes; SY ← ∅; , where l(a) is the range of value domain on a. In this study, we do not deal with missing values. However, it is possible to define µ Ra (x, y) for incomplete information by following the study in [START_REF] Nguyen | On probability of matching in probabilistiy based rough set definitions[END_REF].

SY →X ← ∅; for x ∈ Yi do if µ R Y i (x) ≥ tY ) then SY ← S ∪ {x}; else Z(x) ← j {Xj|µ R X j (x) ≥ tX }; if Z(x) = ∅ then Find Xm ∈ X satisfy µ R Xm (x) = maxj{µ R X j (x)|Xj ∈ Z(x)}; if a pair occurs then choose the first Xm; µX m (x) ← µY i (x); µY i (x) ← 0; SY →X ← SY →X ∪ {x}; end end end S ← i {Xi ∈ X} ∪ SY ∪ SY →X ;
The fuzzy tolerance relation on a set of attributes can be calculated by a combination of the relation on each attribute by Lukasiewicz's t-norm T = max(x 2 + x 1 -1, 0). For those which use OWA fuzzy rough sets, the OWA operators are set as follows:

W min =< w min i >= w min n+1-i = 2 3-i 2 3 -1 for i = 1, 2, 3 and 0 for i = 4, ..., n, W max =< w max i >= w max i = 2 3-i 2 3 -1 for i = 1
, 2, 3 and 0 for i = 4, ..., n, where n is number of instances in datasets. For our approach, we also use OWA fuzzy rough set with above operators and set α = 0.5 to calculate µ R X . That is the same as instances quality measurement in FRIPS-SMOTE-FRBPS in terms of distinguishing instances quality.

Next, selected data are used for training using several classification techniques including k-Nearest Neighbours [START_REF] Cover | Nearest neighbor pattern classification[END_REF], Support Vector Machine [START_REF] Cortes | Support-vector networks[END_REF], naive Bayes [START_REF] Rish | An empirical study of the naive bayes classifier[END_REF] and Decision Trees [START_REF] Strobl | An introduction to recursive partitioning: Rationale, application and characteristics of classification and regression trees, bagging and random forests[END_REF]. It should be noticed that these data classification toolboxes are available within R data mining packages.

To evaluate classification performance in different cases where classifier techniques are used as well as instance selection methods, we use the confusion matrix (Table 2) to calculate the Area Under the ROC Curver(AU C) [START_REF] Huang | Using auc and accuracy in evaluating learning algorithms[END_REF]. AU C provides a single measure for comparing classification on average.

AU C = 1 2
T P T P + F N + T N T N + F P . (15)
For comparison, we try to obtain the best performance for each instance selection algorithm. FRPSS suggests its own method of threshold optimization. For FRIS and MFRISs, instance measurements are quantized into a set of a limited number of elements. We then use this set as threshold candidates in order to find out if there is a significant increase in performance.

In general, from the average AUC shown in Table 3, we can see that after using the modification of rough set instances selection (both MFRIS1 and MRIS2), the classification performance increases noticeably in almost datasets. It seems that the use of multiple thresholds is more effective with the increase of IR value. Although the original fuzzy rough instance selection approach does not deal with imbalanced data, the increase in performance of MFRISs proves that it is possible to use fuzzy rough sets based instance selection with some criteria to overcome the imbalance issue.

Although less consuming time may be an advantage, besides high performance archived, we do not show the time measurements in this version. This is due to the different criteria in choosing thresholds. FRISS, as stated before, has its own method to optimize thresholds. It considers the unique set of all instance measurements as the threshold candidates. This means that the algorithm has to train and test numerous times (nearly as many as the cardinality of training sets). In FRIS and MFRISs, on the other hand, we tried in limited times.

To a comparison of the two algorithms using our approach, MFRIS1 archives better performance for most datasets even though MFRIS2 still improves performance. In the second algorithm, all instances in minority classes are chosen. If such instances do not have sufficient quality, an edited training set from MFRIS2 cannot be better than that from MFRIS1.

When comparing the classifiers, in general, we can see the differences among them in terms of performance evaluation when there is no instance selection method employed (Figure 4, transparency color bars). SVM seems to be sensitive to imbalance issues, while Bayes provides better performance. This means that for a different type of dataset with respect to imbalance issues, a different choice of classifier may be required.

Figure 4 also illustrates the difference in performance increase (dark color bars). SVM does not show a good result for imbalanced data, nevertheless, its performance improves to a greater extent than the rest. In other words, naive Bayes receives less support from instance selection methods. However, the performance of all classifiers generally improves.

Obviously, optimized thresholds are various among classifiers and among datasets. Table 4 shows that SVM needs to remove more instances by using smaller thresholds in comparing with other classifiers. In addition, Table 5 suggests that for some datasets, it is necessary to edit majority classes by discarding many instances (high majority thresholds) while moving some instances to minority classes.

Conclusion and future work

In this paper, we presented a modification of fuzzy-rough instance selection using multiple thresholds. Retaining the advantages of fuzzy-rough sets when selecting instances, it can measure the quality of instances for training phase in order to select valuable instances. It also deals with imbalanced data by using different thresholds between minority and majority classes. Some experiments show that these approaches can improve performance considerably. They can also be competitive with state-of-the-art methods. Furthermore, this notion can potentially applied to modify current rough set based learning methods in order to eliminate the instance selection step.

This paper mainly proves the robustness of using multiple thresholds in fuzzyrough sets based instance selection approach. No algorithm for selecting valuable thresholds is suggested. Therefore, in further research we intend to study a method of optimizing thresholds and will endeavour to discover a better instance quality measurement function.

3 Fig. 1 .

 31 Fig. 1. Histogram to show distribution of approximations and average approximations for a dataset

Fig. 2 .

 2 Fig. 2. MFRIS1 -The first algorithm: Choosing or eliminating instances for both majority and minority classes

Fig. 3 .

 3 Fig. 3. MFRIS2 -The second algorithm: Choosing, removing or relabelling instances for majority classes

Fig. 4 .

 4 Fig. 4. Comparing average of AUC (on all datasets) and increase on average of AUC using MFRISs for different classifiers

Table 1 .

 1 Datasets for experiments

	Datasets	Inst... Feat...	IR Datasets	Inst... Feat...	IR
	pima	768	8 1.87 shuttle-c0-vs-c4	1829	9 13.87
	haberman	305	3 2.77 page-blocks-1-3 vs 4	472	10 15.86
	liver man	439	9 2.78 yeast4	1484	8 28.10
	vehicle1	846	18 2.90 winequality-red-4	1599	11 29.17
	blood transfusion 748	4 3.20 ozone one hr	1848	72 31.42
	life-style-diseases 4819	13 4.04 yeast5	1484	8 32.73
	segment0	2308	19 6.02 ecoli-0-1-3-7 vs 2-6	280	7 39.00
	page-blocks0	5472	10 8.79 yeast6	1484	8 41.40
	yeast-2 vs 4	514	8 9.08 mammography	11183	6 42.01
	vowel0	988	13 9.98 winequality-red-8 vs 6-7	855	11 46.50
	glass2	214	9 11.59 winequality-white-3-9 vs 5 1482	11 58.28
	ozone eight hr	1847	72 13.43 winequality-red-3 vs 5	691	11 68.10

Table 2 .

 2 Confusion matrix for a two-class problem Positive prediction Negative prediction Positive class True Positive (TP) False Negative (FN) Negative class False Positive (FP) True Negative (TN)

Table 3 .

 3 Average AUC archived by the different classifiers for each instance selection method (bold numbers show the highest)

	Datasets	NONE FRIS	FRPSS MFRIS1 MFRIS2
	pima	0.7031	0.7099	0.6967	0.7462 0.7303
	haberman	0.5559	0.5795	0.6292	0.6565 0.6195
	liver man	0.5673	0.5903	0.6665	0.7078 0.6991
	vehicle1	0.6478	0.6545	0.6931	0.7357	0.7387
	blood transfusion	0.5802	0.6084	0.6290	0.6845 0.6765
	life-style-diseases	0.5537	0.5794	0.5729	0.6055 0.5629
	segment0	0.8645	0.8769	0.8638	0.9466 0.9058
	page-blocks0	0.8177	0.8433	0.8650	0.8963 0.8870
	yeast-2 vs 4	0.7892	0.7958	0.7678	0.9102 0.8856
	vowel0	0.8996	0.8443	0.9482	0.9556 0.9498
	glass2	0.5305	0.5489	0.6082	0.7291 0.6918
	ozone eight hr	0.5738	0.5738	0.7325 0.5738	0.5738
	shuttle-c0-vs-c4	0.9947	0.6166	0.9947	0.9987 0.9983
	page-blocks-1-3 vs 4	0.7996	0.8048	0.8059	0.9187 0.9046
	yeast4	0.5710	0.6020	0.8069	0.8416 0.7924
	winequality-red-4	0.5090	0.5094	0.6232	0.6621	0.6785
	ozone one hr	0.5659	0.5659	0.7485 0.5659	0.5659
	yeast5	0.8081	0.8081	0.8736	0.9759 0.9421
	ecoli-0-1-3-7 vs 2-6	0.6487	0.6519	0.5625	0.7766	0.8593
	yeast6	0.6634	0.6879	0.8019	0.8996 0.8587
	mammography	0.7559	0.7735	0.7790	0.8771 0.7587
	winequality-red-8 vs 6-7 0.5085	0.5102	0.5896	0.6661	0.7278
	winequality-white-3-9 vs 5 0.5310	0.5318	0.5715	0.7111	0.7464
	winequality-red-3 vs 5	0.5252	0.5128	0.5385	0.6481	0.8376

Table 4 .

 4 Comparing the thresholds for MFRIS1 that give the highest improvement of each classifier on mammography dataset (the first number is for minority)

	Training sets	kNN	SVM	nBayes	DecisionTree
	fold2 3	0.5090 0.5966 0.5300 0.6483 0.5000 0.5966 0.5090 0.5966
	fold1 3	0.5020 0.5841 0.5220 0.6322 0.4540 0.5841 0.5090 0.6322
	fold1 2	0.5100 0.5911 0.5280 0.6517 0.4040 0.5911 0.5210 0.6517

Table 5 .

 5 Comparing the thresholds for MFRIS2 that the highest improvement of each classifiers on winequality-red-8 vs 6-7 dataset (the first number is for minority)

	Training sets	kNN	SVM	nBayes	DecisionTree
	fold2 3	0.2560 0.9075 0.2560 0.9602 0.1250 0.9360 0.2140 0.9602
	fold1 3	0.2130 0.9085 0.2130 0.8758 0.1600 0.9085 0.1840 0.8758
	fold1 2	0.0000 0.7619 0.1690 0.9575 0.0000 0.9082 0.1440 0.9082

A fuzzy implicator is a function I : [0, 1] × [0, 1] → [0, 1] which satisfies the following properties: I(0, 0) = 1; T (1, a) = a; I is decreasing in the first argument; I is increasing in the second argument.

A t-norm is a function T : [0, 1] × [0, 1] → [0,1] which satisfies the following properties: Commutativity: T (a, b) = T (b, a); Monotonicity: T (a, b) ≤ T (c, d) if a ≤ c and b ≤ d; Associativity: T (a, T (b, c)) = T (T (a, b), c); The number 1 acts as an identity element: T (a, 1) = a

http://www.r-project.org/

http://www.rstudio.com/


Thouraya Merazi-Meksen

Malika Boudraa

Bachir Boudraa

Ultrasonic Data Compression and Analysis during Material Inspection

Keywords: Ultrasounds, Data compression, Signal processing algorithm, Shape recognition

Ultrasonic inspection is nowadays widely used in material inspection to detect and locate internal defects. However, in case of large structures inspection, such as pipelines, vessel reactors etc., the amount of signals to acquire and interpret may be very important, especially in modern techniques that use image representations.

In this article, a proposed algorithm to ultrasonic data selection and processing is described. Images will be replaced by 2 arrays that contain only the data regarding a detected defect. This allows avoiding the storage of useless received ultrasonic data, and simple pattern recognition techniques may be exploited to defect characterization.

Introduction

In ultrasonic examination of materials, Time Of Flight Diffraction (TOFD) technique has demonstrated high accuracy in detection of internal defects presented in materials [START_REF] Silk | The Potential of Scattered or Diffracted Ultrasounds in the Determination of Crack Depth[END_REF]. Based on data digitalization and processing, it improves inspection performances in terms of reliability and sensitivity (better defect detection and identification), especially for crack detection. The Time of Flight of the diffracted ultrasonic echo produced by the extremities of the crack is measured and stored to determine the position of a pair of probes and be displayed as image known as D-SCAN or TOFD image. Since Bolland demonstrated that internal cracks presented in a structure under control involve a set of parabolic shapes in TOFD images [START_REF] Bolland | Imagerie Ultrasonore pour la Detection Automatique de Défauts en Contrôle non Destructif[END_REF], detecting parabolic forms has become a solution to automatically analyze TOFD images for crack detection. Lew Lei et al used a gradient based Hough Transform [START_REF] Lew | Gradient Based Hough Transform for the Detection and Characterisation of Defect during Non Destructive Inspection[END_REF], Maalmi exploited genetic algorithms [4], and Petcher et al. used a matched filtering [START_REF] Petcher | Parabola detection using matched filtering for ultrasound Bscans[END_REF].

When large structures are inspected (reactors, pipelines...) the amount of data to interpret will be considerable, and recording all images, especially when there is no defec,t is useless. The inconvenient of these methods is the processing of all pixels of images, as the most of them are zeros, especially when no defect is presented. This is why recent works aimed to develop automated data interpretation [START_REF] Sambath | Automatic Defect Classification in Ultrasonic NDT Using Artificial Intelligence[END_REF]7], and interest of developing low dimension representations is increasing [START_REF] Zhang | Sparse signal representation and its applications in ultrasonic NDE[END_REF][START_REF] Boβmann | Sparse deconvolution methods for ultrasonic NDT[END_REF].

In this context, we proposed in a precedent work, a method based on neural networks that selects ultrasonic data regarding a detected defect. This allows avoiding the storage and the processing of useless signals [START_REF] Merazi-Meksen | Neural Networks to Select Ultrasonic Data in Non Destructive Testing[END_REF]. The work presented in this paper is another approach to built the sparse matrix structure for ultrasonic data storage in order to record the only ones concerning the defects in one hand, and in another hand, to show that pattern recognition tools such as RANdom Sample Consensus (RANSAC) filter can be applied on the sparse matrix elements to crack recognition, automating thus the detection.

Next section consists in a brief remind of ultrasonic TOFD technique inspection as well as the modeling of the parabolic forms characterizing the crack defect on images. In section 3 the proposed approach for data storage in sparse matrix structure is described as well as the application of the RANSAC filter on its elements. Results and conclusion are given in section 4.

TOFD Inspection

Principle

When ultrasound is incident at linear discontinuity such a crack, diffracted waves take places at its extremities in addition to the normal reflected and refracted waves. This diffracted energy is emitted over a wide angular range and is assumed to originate at the extremities of the flaw (Figure 1). In TOFD inspection, two transducers, are moved step by step according to a straight line, one as emitter and one as receiver (Figure 2). At each position, the received signal is converted in electrical one, digitized, then displayed in grey scale as a row of an image. Each sample of the signal represents a pixel. Usually, zero amplitude light in grey, negative maximum amplitude in black and positive maximum amplitude in white. The figure 3 shows an example. When the defect presented in the materials is a crack, it appears on the image a set of parabolas with a main axis. This zone, called region of interest is relatively small in comparison with the TOFD image. In this work, it is proposed to exploit the sparse matrix structure to register only the pertinent data regarding the detected defect to avoid the record of a huge of zeros (or noise signals) during the control.

Parabola modeling

During a TOFD inspection, the probes are displaced step by step, according to a straight line. When they are at equal distance of the defect (d x = 0), the time of flight defect Transducers Emitter/receivcer Time (samples) Pixels Pixesls t(dx) is minimal. As dx increases, the calculations formulated by Bolland in [START_REF] Bolland | Imagerie Ultrasonore pour la Detection Automatique de Défauts en Contrôle non Destructif[END_REF] conclude that t(dx) is given by Formula (1):

t(d x )) = z 2 c(z 2 +d 2 ) 3/2 (d x ) 2 + 2√z 2 +d 2 c (1)
where c is the ultrasound velocity, z is the depth of the defect and d is half of the distance between the two transducers. Because these three parameters are constants, Formula (1) can be expressed as follows:

t(d x ) = A(d x ) 2 +B (2)
This form is characteristic of a parabola. Figure 5 shows the representation of the time of flight (TOF) according to the displacement d x .

3

Sparse Matrix Construction

Sparse Matrix Structure

The sparse matrix is a compressed manner to store matrix when the major elements are zeros. The principle consists in expressing it by three or two arrays. One array contains the non zeros values, and the two others contain their respective coordinates. For example, the following 5 x 7 matrix M will be represented by the three arrays A, XA, YA two arrays XA and YA are stored. There are two types of sparse matrix: structured and unstructured. In the structured matrix, the nonzero elements form a regular pattern or lie in blocks and a graph may be used for representing it by connecting the elements.

Data Selection

When a matrix contains few non zeros elements, it is termed sparse and may be recorded as a 2D array containing the coordinates on the non zero elements [START_REF] Duff | Direct Method for Sparse Matrix[END_REF]. On ultrasonic images, the zone of interest (defect zone) is very small in comparison of the dimension of the structure under control. It is then suitable to process in the decision phase, only the data regarding the detected defect. The idea consists in avoiding the storage of the signal when it is only constituted by noise, and in another hand to select a unique sample of the reached signal in the defect zone. This sample corresponds to the amplitude determined by a threshold calculated after comparison with the mean and the variance of a signal reached in a zone without defect. Figure 5 shows an example. Thus, calculations will be done on the 120 resulting points, instead of the 150x120 pixels of the TOFD image.

Parabola Detection

The three well known used methods to parameterized form detection are: Least square, the Hough Transform and the RANSAC filter. Least square is the easier algorithm to use if the points belonging to the curve are known. RANSAC is used where many points, (called outliers) do not fit the researched form [START_REF] Fischler | Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography[END_REF]. The points that are nearest the curve are called inliers. The Hough transform is a voting method used when many forms are presented in the image.

In TOFD images, shapes corresponding to defect have parabolic forms. The aim of this work is to detect those parabolas in an automated manner. RANSAC filter is the most suitable algorithm to exploit, as it can detect imperfect parameterized forms (presence of noise, incomplete forms…). This algorithm is simple and need less computational time than the Hough transform. The principle consists in estimating parameters of the shape which can be expressed by a mathematical model with a certain probability. It can thus be applied to detect lines, circles, parabolas on images. The benefit is that imperfect forms can be detected when some data, called outlier, do not fit the model. Outliers can be generated by extreme values of the noise or from incorrect measurements. The points that are the nearest the researched form are called inliers. The figure 7 illustrates an example of a set of data, where the inliers can be fitted to a line. -Calculate a model that uses those 2 points -Determine how many points from the set of all points fi t with a predefi ned tolerance ε .

-Determine the number of inliers Ni over the total number points in the set.

Until Ni > τ inliers outlier

The number of iterations N is calculated by equation 3:

𝑁 = log(1-𝑝) log (1-𝑢) (3)
Where m is the minimum number of points required, u, the probability that any selected data point is an inlier, and p the probability that at least one of the sets of random samples does not include an outlier (usually set to 0.99) .

Results and conclusions

When applying the RANSAC filter on the sparse matrix represented in figure 6(b), the calculated parabola is drawn to illustrate the result (figure 8.). The three prominent points are examples of randomly selected samples to parabola calculation.

Fig. 8. The optimal calculated parabola is drawn in continuous line. The three prominent points are an example of a set of 3 points randomly selected from the selected samples of data constituting the sparse matrix, built instead of TOFD image on figure 6(a). In precedent works [START_REF] Bolland | Imagerie Ultrasonore pour la Detection Automatique de Défauts en Contrôle non Destructif[END_REF][START_REF] Lew | Gradient Based Hough Transform for the Detection and Characterisation of Defect during Non Destructive Inspection[END_REF][4][START_REF] Petcher | Parabola detection using matched filtering for ultrasound Bscans[END_REF] parabolas were respectively detected by calculating the gradient, matched filters and genetic algorithms. All pixels have to be considered in the processing what involves time consuming and useless data storage.

This work shows the feasibility of avoiding TOFD image formation, replacing it by a sparse matrix on which RANSAC filter may be applied. First, a pre-processing selects on every row of the region of interest, a single pixel, sufficient to determine the time of flight of the signal echo at corresponding probes position. Recording only those two coordinates (Time of Flight and probe position) during the probe displacement led to a matrix formation. TOFD image is thus reduced to a sparse matrix and the points to consider considerably reduced. In the decision phase, RANSAC filter had been exploited in order to recognize the parabolic form and determining its summits coordinates to crack location.

Fig. 1 .

 1 Fig. 1 .Ultrasonics waves.

1 :Fig. 2 .

 12 Fig. 2. Structure under control for defect detection.

Fig. 3 .

 3 Fig. 3. (a) Ultrasonic signal, (b) TOFD image.

Fig. 4 .

 4 Fig. 4. Theoretical representation of the relationship between the time of flight t and the displacement dx of the transducers.

  There are two types of sparse matrix: Many zeros elements or many zeros lines. When the only non zeros positions values are needed (as in this studied case), the only A = [1 3 9 2 4] (non zero values) XA = [2 2 3 4 4] (corresponding rows) YA = [2 3 2 3 4] (corresponding columns)

Fig. 5 .

 5 Fig. 5. (a) Ultrasonic reached signal, (b) Selected sample corresponding to the Time of Flight.

Figure 6 (Fig. 6 .

 66 Figure 6 (b) shows the resulted sparse matrix that replaces the image of figure 6(a).

Fig. 7 .

 7 Fig. 7. Example of points fitting a line.

Figure 9 (Fig. 9 .

 99 Figure 9 (a) represents another example of a TOFD image presenting three detected cracks. In (b), stars represent sparse matrix elements resulting and the parabolas calculated using the RANSAC algorithm.


Andrés Iglesias
email: iglesias@unican.es

Akemi Gálvez

Marta Collantes

Global-Support Rational Curve Method for Data Approximation with Bat Algorithm

Keywords: data approximation, rational curve, bat algorithm, metaheuristic technique, Bézier curve

The problem of obtaining an approximating curve from a given set of data points appears recurrently in several applied and industrial domains, such as CAD/CAM, computer graphics and animation, medicine, and many others. Although polynomial blending functions are usually applied to tackle this issue, some shapes cannot yet be adequately approximated by using the polynomial scheme. In this paper we address this limitation by applying rational global-support blending functions, particularly rational Bézier curves. Our method is based on a natureinspired meta-heuristic called bat algorithm, which has been recently introduced to solve difficult optimization problems. To check the performance of our approach, it has been applied to some illustrative examples of 2D and 3D curves. Our results show that the method performs very well, being able to yield a satisfactory approximating curve with a high degree of accuracy.

Introduction

Motivation

Obtaining an approximating curve from a given set of data points is a very important and recurrent problem in many applied and industrial domains. It appears, for instance, in computer-aided design and manufacturing (CAD/CAM), a field where data points are usually obtained from real measurements of an existing geometric entity, as it typically happens in the construction of car bodies, ship hulls, airplane fuselage, and other free-form objects. It is also a common problem in the shoes industry, in archeology (reconstruction of archeological assets), in medicine (computer tomography (CT), magnetic resonance imaging), computer graphics and animation, virtual and augmented reality, and in many other fields.

Data points in these settings are acquired through many different technologies, such as 3D laser scanning, touch scanners, coordinate measuring machines, CT scanners, or convergent photogrammetry, to mention just a few. A common factor of all these techniques is that a huge number of data (often in the range of hundreds of thousands, and even millions) is usually obtained. Besides, in most cases no geometric or topological information is available beyond the data points. Since dealing with this amount of data becomes impractical, a primary goal in the field is to convert the real data from a physical object into a fully usable digital model, a process commonly called reverse engineering. This allows significant savings in terms of storage capacity and processing and manufacturing time. Furthermore, the digital models are easier and cheaper to modify than their real counterparts and are usually available anytime and anywhere.

In real-world problems, data points are usually affected by measurement noise, irregular sampling, and other artifacts [START_REF] Barnhill | Geometric Processing for Design and Manufacturing[END_REF][START_REF] Patrikalakis | Shape Interrogation for Computer Aided Design and Manufacturing[END_REF][START_REF] Pottmann | Industrial geometry: recent advances and applications in CAD[END_REF]. Consequently, a good fitting of data should be generally based on approximation schemes rather than on interpolation. In this case, the approximating curve is not required to pass through all input data points, but just near to them, according to some prescribed distance criteria.

Several approximating families of functions have been applied to this problem. Among them, the free-form parametric curves such as Bézier, B-spline and NURBS, are widely applied in many industrial settings due to their great flexibility and the fact that they can represent smooth shapes with only a few parameters [START_REF] Barnhill | Geometric Processing for Design and Manufacturing[END_REF][START_REF] Jing | Fitting B-spline curves by least squares support vector machines[END_REF][START_REF] Li | Adaptive knot placement in B-spline curve approximation[END_REF][START_REF] Park | An error-bounded approximate method for representing planar curves in B-splines[END_REF][START_REF] Park | B-spline curve fitting based on adaptive curve refinement using dominant points[END_REF]. In general, the approximating curves can be classified as global-support and local-support. By global-support curves we mean curves mathematically expressed as a combination of basis functions whose support is the whole domain of the problem. As a consequence, these curves exhibit a global control, in the sense that any modification of the shape of the curve in a particular location is propagated throughout the whole curve. This is in clear contrast to the local-support approaches, which provide local control of the shape of the curve [START_REF] Loucera | Simulated annealing algorithm for Bézier curve approximation[END_REF][START_REF] Pottmann | Industrial geometry: recent advances and applications in CAD[END_REF] and have become prevalent in CAD/CAM and computer graphics. In this work we focus on the global-support approach. Some previous papers addressed this problem by using Bézier curves [START_REF] Gálvez | Cuckoo search with Lévy flights for weighted Bayesian energy functional optimization in global-support curve data fitting[END_REF][START_REF] Loucera | Simulated annealing algorithm for Bézier curve approximation[END_REF], which are given by a linear combination of polynomial basis functions (the Bernstein polynomials). Although they obtained good results for a number of shapes, this polynomial approach is still limited, as it cannot adequately describe some particular shapes (such as the conics). As a consequence, there is still a need for more powerful blending functions.

An interesting extension in this regard is given by the rational basis functions, which are mathematically described as the quotient of two polynomials. A remarkable advantage of this rational scheme is that the conics can be canonically described as rational functions. In this paper, we take advantage of this valuable feature to solve the curve approximation problem by using rational Bézier curves. Unfortunately, this rational approach becomes more difficult than the polynomial one, since new parameters are now introduced into the problem. Consequently, we are confronted with the challenge of obtaining optimal values for many (qualitatively different) parameters, namely, data parameters, poles, and weights. This leads to a difficult over-determined multivariate continuous nonlinear optimization problem.

Aims and Structure of the Paper

In this paper, we solve this optimization problem by applying a recently proposed nature-inspired meta-heuristic technique called bat algorithm, which is receiving increasing attention from the scientific community during the last few years due to its good performance in solving difficult continuous optimization problems [START_REF] Yang | A new metaheuristic bat-inspired algorithm[END_REF]. This algorithm is based on the echolocation behavior of microbats with varying pulse rates of emission and loudness (see Section 3 for details). We would like to remark that, in spite of its valuable features for continuous optimization, to the best of authors' knowledge the bat algorithm has never been applied in the context of data fitting for geometric modeling or computer graphics. The present paper aims to fill this gap. In this work, the bat algorithm is applied to obtain a very accurate fitting curve to a given set of data points by using rational global-support blending functions. To check the performance of our approach, it has been applied to three simple yet illustrative examples of two-dimensional and three-dimensional shapes.

The structure of this paper is as follows: previous work in the subject of data fitting with free-form parametric curves is briefly reported in Section 2. In Section 3 we provide a gentle overview about the bat algorithm and its main features and advantages. The problem of curve approximation with rational Bézier curves along with our proposed approach to solve it are described in Section 4. The section also discusses the important issue of parameter tuning. Then, three illustrative examples of its application are reported in Section 5. Our experimental results show that the presented method performs very well, being able to replicate the underlying shape of data accurately. The paper closes with the main conclusions of this contribution and our plans for future work in the field.

Previous Work

The problem of data approximation with free-form parametric curves has been the subject of research for many years. First approaches in the field were mostly based on numerical procedures [START_REF] Dierckx | Curve and Surface Fitting with Splines[END_REF][START_REF] Farin | Curves and surfaces for CAGD[END_REF][START_REF] Powell | Curve fitting by splines in one variable[END_REF]. However, it has been shown that traditional mathematical optimization techniques fail to solve the problem in its generality. Consequently, there has been a great interest to explore other possible approaches to this problem. Some recent approaches in this line use error bounds [START_REF] Park | An error-bounded approximate method for representing planar curves in B-splines[END_REF], curvature-based squared distance minimization [START_REF] Wang | Fitting B-spline curves to point clouds by curvature-based squared distance minimization[END_REF], or dominant points [START_REF] Park | B-spline curve fitting based on adaptive curve refinement using dominant points[END_REF]. In general, they perform well but require some particular constraints (such as high differentiability, closed curves, noiseless data) which are not so commonly met in real-world applications.

On the other hand, interesting research carried out during the last two decades has shown that the application of artificial intelligence and soft computing techniques can achieve remarkable results for this problem [START_REF] Barhak | Parameterization and reconstruction from 3D scattered points based on neural network and PDE techniques[END_REF][START_REF] Gu | Neural network approach to the reconstruction of free-form surfaces for reverse engineering[END_REF][START_REF] Hoffmann | Numerical control of Kohonen neural network for scattered data approximation[END_REF]. Most of these methods rely on some kind of neural networks, such as standard neural networks [START_REF] Gu | Neural network approach to the reconstruction of free-form surfaces for reverse engineering[END_REF], or Kohonen's SOM (Self-Organizing Maps) nets [START_REF] Hoffmann | Numerical control of Kohonen neural network for scattered data approximation[END_REF]. In some cases, this neural approach is combined with partial differential equations [START_REF] Barhak | Parameterization and reconstruction from 3D scattered points based on neural network and PDE techniques[END_REF] or other approaches [START_REF] Iglesias | Hybrid functional-neural approach for surface reconstruction[END_REF]. The generalization of these methods to functional networks is also analyzed in [START_REF] Echevarría | Extending neural networks for B-spline surface reconstruction[END_REF][START_REF] Iglesias | Functional networks for B-spline surface reconstruction[END_REF]. The application of support vector machines to solve the least-squares B-spline curve fitting problem is reported in [START_REF] Jing | Fitting B-spline curves by least squares support vector machines[END_REF].

Other approaches are based on the application of nature-inspired metaheuristic techniques, which have been intensively applied to solve difficult optimization problems that cannot be tackled through traditional optimization algorithms. Genetic algorithms have been applied to this problem in both the discrete version [START_REF] Sarfraz | Capturing outline of fonts using genetic algorithms and splines[END_REF] and the continuous version [START_REF] Gálvez | Iterative two-step genetic-algorithm method for efficient polynomial B-spline surface reconstruction[END_REF][START_REF] Yoshimoto | Data fitting with a spline using a realcoded algorithm[END_REF]. Other metaheuristic approaches applied to this problem include the use of the popular particle swarm optimization technique [START_REF] Gálvez | Efficient particle swarm optimization approach for data fitting with free knot B-splines[END_REF][START_REF] Gálvez | Particle swarm optimization for non-uniform rational Bspline surface reconstruction from clouds of 3D data points[END_REF], artificial immune systems [START_REF] Gálvez | Immunological-based approach for accurate fitting of 3D noisy data points with Bézier surfaces[END_REF][START_REF] Gálvez | Elitist clonal selection algorithm for optimal choice of free knots in B-spline data fitting[END_REF], firefly algorithm [START_REF] Gálvez | Firefly algorithm for polynomial Bézier surface parameterization[END_REF][START_REF] Gálvez | From nonlinear optimization to convex optimization through firefly algorithm and indirect approach with applications to CAD/CAM[END_REF][START_REF] Gálvez | Firefly algorithm for explicit B-Spline curve fitting to data points[END_REF], cuckoo search [START_REF] Gálvez | Cuckoo search with Lévy flights for weighted Bayesian energy functional optimization in global-support curve data fitting[END_REF], simulated annealing [START_REF] Loucera | Simulated annealing algorithm for Bézier curve approximation[END_REF], estimation of distribution algorithms [START_REF] Zhao | Adaptive knot adjustment using a GMM-based continuous optimization algorithm in B-spline curve approximation[END_REF], memetic algorithms [START_REF] Gálvez | New memetic self-adaptive firefly algorithm for continuous optimization[END_REF], and hybrid techniques [START_REF] Gálvez | A new iterative mutually-coupled hybrid GA-PSO approach for curve fitting in manufacturing[END_REF][START_REF] Sarfraz | Capturing outline of fonts using genetic algorithms and splines[END_REF]. These methods yield good results for shapes that can be properly described in polynomial terms, but more complicated shapes are still elusive. Furthermore, it has been shown that even some simple shapes such as the conics cannot be adequately described in terms of polynomial functions. In this paper we overcome this limitation by considering rational basis functions, as discussed in section 4.

The Bat Algorithm

Basic Principles

The bat algorithm is a bio-inspired population-based meta-heuristic algorithm originally proposed by Xin-She Yang in 2010 to solve difficult optimization problems [START_REF] Yang | A new metaheuristic bat-inspired algorithm[END_REF][START_REF] Yang | Bat algorithm: a novel approach for global engineering optimization[END_REF]. The algorithm is based on the echolocation behavior of bats. The author focused particularly on microbats, as they use a type of sonar called echolocation, with varying pulse rates of emission and loudness, to detect prey, avoid obstacles, and locate their roosting crevices in the dark.

Despite the short time since its appearance, the bat algorithm has already been applied to several engineering and industrial problems. Simultaneously, some modifications and improvements on the original version have been developed, such as the muti-objective bat algorithm [START_REF] Yang | Bat algorithm for multiobjective optimization[END_REF], directed artificial bat algorithm [START_REF] Rekaby | Directed artificial bat algorithm (DABA): a new bio-inspired algorithm[END_REF], binary bat algorithm [START_REF] Mirjalili | Binary bat algorithm[END_REF], and others. The interested reader is referred to the general paper in [START_REF] Yang | Bat algorithm: literature review and applications[END_REF] for a comprehensive, updated review and taxonomic classification of the bat algorithm and all its variants and applications.

In this paper we consider the standard bat algorithm, as described in the original paper in [START_REF] Yang | A new metaheuristic bat-inspired algorithm[END_REF]. According to that source, the idealization of the echolocation of microbats can be summarized as follows:

1. Bats use echolocation to sense distance and distinguish between food, prey, and background barriers. 2. Each virtual bat flies randomly with a velocity v i at position (solution)

x i with a fixed frequency f min , varying wavelength λ and loudness A 0 to search for prey. As it searches and finds its prey, it changes wavelength (or frequency) of their emitted pulses and adjust the rate of pulse emission r, depending on the proximity of the target. 3. It is assumed that the loudness will vary from an (initially large and positive) value A 0 to a minimum constant value A min .

In order to apply the bat algorithm for optimization problems more efficiently, some additional assumptions are strongly advisable. In general, we assume that the frequency f evolves on a bounded interval [f min , f max ]. This means that the wavelength λ is also bounded, because f and λ are related to each other by the fact that the product λ.f is constant. For practical reasons, it is also convenient that the largest wavelength is chosen such that it is comparable to the size of the domain of interest (the search space, for optimization problems). For simplicity, we can assume that f min = 0, so f ∈ [0, f max ]. The rate of pulse can simply be in the range r ∈ [0, 1], where 0 means no pulses at all, and 1 means the maximum rate of pulse emission. With these idealized rules indicated above, the basic pseudo-code of the bat algorithm is shown in Algorithm 1. It is described in next section.

The Algorithm

Basically, the algorithm considers an initial population of P individuals (bats).

Each bat, representing a potential solution of the optimization problem, has a location x i and velocity v i . The algorithm initializes these variables with random values within the search space. Then, the pulse frequency, pulse rate, and loudness are computed for each individual bat (lines 2-6). Then, the swarm evolves in a discrete way over generations (line 7), like time instances (line 21) until the maximum number of generations, G max , is reached (line 22). For each generation g and each bat (line 8), new frequency, location and velocity are computed (lines 9-10) according to the following evolution equations:

f g i = f g min + β(f g max -f g min ) (1)
v g i = v g-1 i + [x g i -x * ] f g i (2)
x g i = x g-1 i + v g i (3)
where β ∈ [0, 1] follows the random uniform distribution, and x * represents the current global best location (solution), which is obtained through evaluation of the objective function at all bats and ranking of their fitness values. The superscript (.) g is used to denote the current generation g. It is worthwhile to remark a certain similarity of the evolutions equations ( 2)-(3) to those of the velocity and position in standard particle swarm optimization (PSO), as f i controls the pace and range of the movement of the particles of the swarm. From this viewpoint, the bat algorithm is a balanced combination of the standard PSO and the local search modulated by the loudness and pulse rate.

The current global best solution and a local solution around it are probabilistically selected according to a given criterion (lines 11-14). Then, search is Initialize pulse rates ri and loudness Ai 6: end for 7: while g ≤ Gmax do 8:

for i = 1 to P do 9:

Generate new solutions by adjusting frequency, 10:

and updating locations and velocities

16:

if ψi < Ai and φ(x i ) < φ(x * ) then 17:

Accept new solutions 18:

Increase ri and decrease Ai //eqns. ( 5)-( 6)

19: end if 20:

end for 21:

g ← g + 1 22: end while 23: Rank the bats and find current best x * 24: return x * Algorithm 1: Bat Algorithm intensified by a local random walk (line 15). For this local search, once a solution is selected among the current best solutions, it is perturbed locally through a random walk of the form:

x new = x old + A g (4)
where is a random number with uniform distribution on the interval [-1, 1] and A g =< A g i >, is the average loudness of all the bats at generation g. If the new solution achieved is better than the previous best one, it is probabilistically accepted depending on the value of the loudness. In that case, the algorithm increases the pulse rate and decreases the loudness (lines [START_REF] Gálvez | Iterative two-step genetic-algorithm method for efficient polynomial B-spline surface reconstruction[END_REF][START_REF] Gu | Neural network approach to the reconstruction of free-form surfaces for reverse engineering[END_REF][START_REF] Hoffmann | Numerical control of Kohonen neural network for scattered data approximation[END_REF][START_REF] Iglesias | Functional networks for B-spline surface reconstruction[END_REF]. This process is repeated for the given number of generations. In general, the loudness decreases once a bat finds its prey (in our analogy, once a new best solution is found), while the rate of pulse emission decreases. For simplicity, the following values are commonly used: A 0 = 1 and A min = 0, assuming that this latter value means that a bat has found the prey and temporarily stop emitting any sound. However, any other value within its range can also be feasible (see our discussion about parameter tuning in Section 4.2 for details). The evolution rules for loudness and pulse rate are as follows:

A g+1 i = αA g i (5)
r g+1 i = r 0 i [1 -exp(-γg)] (6)
where α and γ are constants. Note that for any 0 < α < 1 and any γ > 0 we have:

A g i → 0, r g i → r 0 i , as g → ∞ (7)
In general, each bat should have different values for loudness and pulse emission rate, which can be computationally achieved by randomization. To this aim, we can take an initial loudness A 0 i ∈ (0, 2) while the initial emission rate r 0 i can be any value in the interval [0, 1]. Loudness and emission rates will be updated only if the new solutions are improved, an indication that the bats are moving towards the optimal solution. As a result, the bat algorithm applies a parameter tuning technique to control the dynamic behavior of a swarm of bats. Similarly, the balance between exploration and exploitation can be controlled by tuning algorithm-dependent parameters.

The Proposed Method

Problem To Be Solved

We assume that the reader is familiar with the main concepts of free-form parametric curves [START_REF] Farin | Curves and surfaces for CAGD[END_REF]. Mathematically, a free-form rational Bézier curve Φ(τ ) of degree η is defined as:

Φ(τ ) = η j=0 ω j Λ j φ η j (τ ) η j=0 ω j φ η j (τ ) (8)
where Λ j are vector coefficients called the poles, ω j are their scalar weights, φ η j (τ ) are the Bernstein polynomials of index j and degree η, given by:

φ η j (τ ) = η j τ j (1 -τ ) η-j
and τ is the curve parameter, defined on the finite interval [0, 1]. By convention, 0! = 1. Note that in this paper vectors are denoted in bold. Suppose now that we are given a set of data points {∆ i } i=1,...,κ in R ν (usually ν = 2 or ν = 3). Our goal is to obtain the rational Bézier curve Φ(τ ) performing discrete approximation of the data points {∆ i } i . To do so, we have to compute all parameters (i.e. poles Λ j , weights ω j , and parameters τ i associated with data points ∆ i , for i = 1, . . . , κ, j = 0, . . . , η) of the approximating curve Φ(τ ) by minimizing the least-squares error, Υ , defined as the sum of squares of the residuals:

Υ = minimize {τ i } i {Λ j } j {ωj }j        κ i=1        ∆ i - η j=0 ω j Λ j φ η j (τ i ) η j=0 ω j φ η j (τ i )        2        . ( 9
)
Now, taking:

ϕ η j (τ ) = ω j φ η j (τ ) η k=0 ω k φ η k (τ ) (10)
Eq. ( 9) becomes:

Υ = minimize {τ i } i {Λ j } j {ωj }j    κ i=1   ∆ i - η j=0 Λ j ϕ η j (τ )   2    , (11)
which can be rewritten in matrix form as:

Ω.Λ = Ξ (12)
called the normal equation, where:

Ω = [Ω i,j ] =   κ k=1 ϕ η i (τ k )ϕ η j (τ k ) i,j   , Ξ = [Ξ j ] =   κ k=1 ∆ k ϕ η j (τ k ) j 
 , Λ = (Λ 0 , . . . , Λ η ) T , for i, j = 0, . . . , η, and

(.) T means the transposition of a vector or a matrix. In general, κ >> η meaning that the system ( 12) is over-determined. If values are assigned to the τ i , our problem can be solved as a classical linear least-squares minimization, with the coefficients {Λ i } i=0,...,η as unknowns. This problem can readily be solved by standard numerical techniques. On the contrary, if the values of τ i are treated as unknowns, the problem becomes much more difficult. Indeed, since the polynomial blending functions φ η j (τ ) are nonlinear in τ and so are the rational blending functions ϕ η j (τ ), the least-squares minimization of the errors is a nonlinear continuous optimization problem. Note also that in many practical cases the number of data points can be extremely large, meaning that we have to deal with a large number of unknowns. In other words, we are also confronted with a high-dimensional problem. It is also a multimodal problem, since there might be arguably more than one set of parameter values leading to the optimal solution.

In conclusion, it is clear that the interplay among all sets of unknowns (data parameters, poles, and weights) leads to a very difficult over-determined, multimodal, multivariate, continuous, nonlinear optimization problem. In this work, we are interested to solve this general problem. This means that we make no assumption about the values of the free parameters; instead, we compute them at full extent.

Proposed Method and Parameter Tuning

Our approach to solve this general problem consists of applying the bat algorithm described above to determine suitable parameter values for the least-squares minimization of functional Υ according to [START_REF] Gálvez | Firefly algorithm for polynomial Bézier surface parameterization[END_REF]. To this aim, each bat, representing a potential solution, corresponds to a parametric vector S j of length κ + η + 1 given by S j = [T j ; W j ], where T j = (τ j 1 , τ j 2 , . . . , τ j κ ) ∈ [0, 1] κ with the {τ j i } i=1,...,κ strictly increasing parameters, and W j = (w j 0 , . . . , w j η ) for j = 1, . . . , P. These parametric vectors are initialized with random values. Also, coordinates of T j are normalized and then sorted. Regarding W j , in this paper we consider values for the weights within the range (0, 100), as values larger than 100 do not modify the shape of the curve noticeably. Application of the bat algorithm described above yields new positions and velocities of the bats representing the potential solutions of our optimization problem.

A critical issue when working with meta-heuristic techniques is the parameter tuning. It is well-known that the performance of these methods is strongly dependent on the choice of suitable values for their parameters. Moreover, such values are problem-dependent, making it hard to determine good values in advance. Therefore, although there are some papers describing suitable values for some problems, our choice must be necessarily empirical. To this purpose, we carried out numerous computer simulations for different parameter values. The different parameters used in our method are arranged in rows in Table 1. For each parameter, the table shows (in columns) its symbol, meaning, range of values, and the parameter value chosen in this paper.

The most important parameters in the bat algorithm are:

population size: in general, increasing the number of individuals (bats) decreases the number of required iterations, but it also increases the number of function evaluations. Therefore, a trade-off between both situations must be achieved for better performance. In this work, we tested populations ranging from 10 to 200 bats and found that while very low values require too many iterations, values in the range 100 -200 perform similarly, so we set this value to 100 individuals. maximum number of iterations: we tested our method for values of this parameter in the range 100 -3000 and found that the method converged in less than 1500 iterations in all our executions. We finally set this parameter in 2000 iterations, but to prevent wasting computation time without any improvement, also set an additional termination criterion: the method stops if no further improvement of the solution is reached after 20 consecutive iterations, even although the total number of iterations is less than 2000. With this additional criterion, the computation times improved significantly without penalizing the quality of the final solution.

initial and minimum loudness and parameter α: they are set to 0.5, 0, and 0.6, respectively. However, from our computer experiments we noticed that our results do not change significantly for values of the initial loudness in the whole range (0, 2), meaning that this parameter is very robust against variations on that interval. Empirically, we found α = 0.6 to be a suitable value for this problem. initial pulse rate and parameter γ: of these two parameters, the initial pulse rate is the most relevant. In fact, parameter γ only affects the very early iterations. We set the initial pulse rate to 0.5, meaning that the selection has an equal probability of change in the long term.

With this choice of parameters, the bat algorithm is run iteratively for a given number of generations or until the convergence of the minimization of the error is eventually achieved. Final positions and velocities of the bats are computed and ranked according to our fitness function Υ . The position of the global best is taken as the final solution of our minimization problem.

Experimental Results

The method described in previous section has been applied to several examples chosen by the authors. Unfortunately, the field still lacks a standardized benchmark for further analysis and comparative purposes. We think, however, that the examples reported here will be useful to determine the good applicability of our method to this problem. To keep the paper in manageable size, in this section we describe only three of them, corresponding to 2D and 3D curves as described in Table 2. In this table, the three examples are reported in rows. For each example, the table shows (in columns) the number of data points used in our experiments along with some other interesting features: if the curve is open or closed, planar or 3D, or whether it includes any non-differentiable point (such as cusps or discontinuities). It also reports the number of degrees of freedom (DOFs) of the optimization problem (i.e. the number of variables to be minimized). These examples have been primarily chosen to reflect the diversity of situations our method can be applied to. First example corresponds to a planar closed curve called hypocycloid, which has several cusps; second example shows a spiral, a planar open curve with a smooth yet challenging shape; and last example corresponds to a helix, a 3D open curve with several DOFs.

Figure 1 shows (from top to bottom) the results of our experiments for the three examples in our benchmark. In all cases, the original data points are displayed as red empty circles while the reconstructed points are displayed as blue + symbols. Note the very good matching between both sets of points for the three examples. This good visual behavior is confirmed by our numerical results, reported in Table 3. The three examples are arranged in rows. For each example, the table reports (in columns):

the average and the best error of the coordinates, according to the equation:

Av µ = κ i=1 |δ µ j -δµ j |
, where δ µ j and δµ j represent respectively the µ coordinate of the input data and reconstructed data, with µ = x, y for 2D curves (µ = x, y, z for 3D curves); the average and best error of the functional Υ according to Eq. ( 9); the average and best error of the root-mean square error, given by:

RM SE = Υ κ .
The average error has been obtained as the mean value from 50 independent executions of the algorithm, while the best value corresponds to the results of the best execution from the 50 runs. As the reader can see, average errors are of order 10 -4 (average) to 10 -5 (best) for the Υ error, with RMSE of order 10 -4 (average) to 10 -5 (best) for the three examples in our benchmark. Coordinate errors show more noticeable variations for the third example due to the scale factor of the vertical component. From the data in Table 3 we can conclude that the method performs very well, being able to replicate the original shape with high accuracy for all instances in our benchmark.

Conclusions and Future Work

This paper introduces a new method for discrete approximation of data points with rational Bézier curves. Given a set of data points, the method computes all relevant parameters (poles, weights, and data parameters) of the rational Bézier fitting curve as the solution of a challenging over-determined nonlinear multimodal multivariate continuous optimization problem. The method is based on a nature-inspired meta-heuristic called bat algorithm, recently introduced to solve difficult optimization problems. To check the performance of our approach, it has been applied to some illustrative examples of 2D and 3D curves. Our results show that the method performs very well, being able to yield a satisfactory approximating curve with a high degree of accuracy. This approach generalizes previous methods for data fitting based on polynomial basis functions to rational blending functions, thus expanding the potential range of applications to include more difficult shapes.

Regarding the implementation issues, all computations in this paper have been performed on a 2.6 GHz. Intel Core i7 processor with 8 GB. of RAM. The source code has been implemented by the authors in the native programming language of the popular scientific program Matlab, version 2013b. The systems of equations have been computed by using specialized Matlab commands for Gaussian elimination with partial pivoting and singular value decomposition (SVD) for squared and non-squared systems, respectively.

Our future work includes the extension of this method to the case of surfaces. We are also interested to analyze its application to some industrial processes and other interesting real-world problems in different fields. Finally, a theoretical analysis about the convergence of this method and a comparative analysis with other alternative approaches on a standardized benchmark (when available) are also part of our future goals. and the University of Cantabria (Santander, Spain). The authors are particularly grateful to the Department of Information Science of Toho University for all the facilities given to carry out this work.

  Require: (Initial Parameters) Population size: P Maximum number of generations: Gmax Loudness: A Pulse rate: r Maximum frequency: fmax Dimension of the problem: d Objective function: φ(x), with x = (x1, . . . , x d ) T Random vectors: Θ = (θ1, . . . , θP ), Ψ = (ψ1, . . . , ψP ) with θ k , ψ k ∈ U (0, 1) 1: g ← 0 //g: generation index 2: for i = 1 to P do 3: Initialize the location and velocity xi and vi //Initialization phase 4: Define pulse frequency fi at xi 5:

Fig. 1 .

 1 Fig. 1. Examples of application of the bat algorithm for data approximation with rational Bézier curves: (top) hypocycloid; (middle) spiral; (bottom) helix. Original points are displayed as red empty circles and the reconstructed points as blue + symbols.

Table 1 .

 1 Parameters and values used in our method.

			Range of Chosen
	Symbol	Meaning	Values	Value
	P	population size	10 -200	100
	g	maximum number 100 -3000 2000 of generations
	A 0	initial loudness	(0, 2)	0.5
	Amin minimum loudness	[0, 1]	0
	r 0	initial pulse rate	[0, 1]	0.5
	fmax maximum frequency [0, 10]	2
	α	multiplicative factor	(0, 1)	0.6
	γ	exponential factor	[0, 1]	0.4

Table 2 .

 2 Benchmark used in this paper along with the main features of each example.

		# data Open/closed Non-differentiable 2D/3D DOFs (degrees
		curve	curve	points	curve of freedom)
	Hypocycloid 200	closed		2D	124
	Spiral	100	open	×	2D	136
	Helix	100	open	×	3D	242

Table 3 .

 3 Fitting errors for the examples used in this paper (arranged in rows): coordinate error, Υ error and RM SE for the mean and best results from 50 executions (in columns).

	Curve	Avµ(mean)	Avµ(best)	Υ (mean)	Υ (best) RMSE (mean) RMSE (best)
	Hypocycloid	x : 6.390788e-4 x : 3.294507e-5 3.475935e-4 3.075698e-5 1.318322e-3 3.921542e-4 y : 7.138176e-4 y : 4.051993e-5
	Spiral	x : 6.998090e-3 x : 2.532214e-4 8.589174e-4 1.867033e-5 2.930729e-3 4.320918e-4 y : 7.189277e-3 y : 2.985567e-4
		x : 1.050018e-2 x : 3.277831e-3
	Helix	y : 1.102567e-2 y : 2.689946e-3	3.261139e-4 1.860439e-5 1.217693e-3 4.313281e-4
		z : 9.711053e-4 y : 3.091722e-4

Acknowledgements. This research has been kindly supported by the Computer Science National Program of the Spanish Ministry of Economy and Competitiveness, Project Ref. #TIN2012-30768, Toho University (Funabashi, Japan),