Lina Hong

Rong Chen

Statistical Fault Localization with Reduced Program Runs

Keywords: Software Fault Localization, Path Redundancy, Statistical Method, Clustering

A typical approach to software fault location is to pinpoint buggy statements by comparing the failing program runs with some successful runs. Most of the research works in this line require a large amount of failing runs and successful runs. Those required execution data inevitably contain a large number of redundant or noisy execution paths, and thus leads to a lower efficiency and accuracy of pinpointing. In this paper, we present an improved fault localization method by statistical analysis of difference between reduced program runs. To do so, we first use a clustering method to eliminate the redundancy in execution paths, next calculate the statistics of difference between the reduced failing runs and successful runs, and then rank the buggy statements in a generated bug report. The experimental results show that our algorithm works many times faster than Wang's, and performs better than competitors in terms of accuracy.

Introduction

A typical thinking in fault localization is to compare successful runs and failing runs [START_REF] Renieris | Fault localization with nearest neighbor queries[END_REF][START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF][START_REF] Wang | Automated path generation for software fault localization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF][START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF]. There are different ways of comparison, and can be divided into distance measures-based methods [START_REF] Renieris | Fault localization with nearest neighbor queries[END_REF][START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF][START_REF] Wang | Automated path generation for software fault localization[END_REF] and characteristic statistics-based methods [START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF][START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF]. But both methods require a lot of failing runs and successful runs. Our analysis found that many characteristics of the runs are same. Figure 1.1 shows a program example with a bug in the assignment statement on line 11, which should be "max = x;" In order to locate this buggy statement by using Wang [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF] method, we run manual design of 12 test cases, and obtain 8 successful runs and 4 failing runs, in which there exists exactly same paths, such as max (4, 2, 3) and max [START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF]. To be fair, for sorting methods students write, we execute randomly generated test cases, and statistics showed that the proportion of redundant paths is so much as 20% ~50%.

Redundancy leads to many problems: [START_REF] Jones | Visualization of test information to assist fault localization[END_REF] Redundancy contributes nothing to the fault location calculations. If two failing runs are with the same path, it means that Fig. 1.A small program they contributions same to fault localization, so there is no need to separately calculate the difference between each of them and successful runs, the same applies two successful runs with two paths the same. (2) Computational efficiency of the system is reduced. Because the system may spend a lot of time doing pointless things.

(3) The calculation of too much redundant data may lead to biased results. Thus removing the path redundancy is necessary.

Measures-based method returns executing differences as bug report, however the statements in differences are not ordered by their importance in literature [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF]. Techniques based on characteristics statistics rank the statements or predicates in order of suspicious, but it is so many statements and predicates, even including those are not predictive of anything which is usually a majority, statistical work on them are time-consuming and effortless.

To overcome these shortcomings, in this paper, we first propose a clustering method to eliminate the path redundancy, and then through statistical analysis of differences between classes of runs to get and rank suspicious statements. Experimental comparison with Wang's algorithm shows that our algorithm makes a greatly improved efficiency and accuracy of locating. In addition, our algorithm can not only diagnose the reasons for each failing run, it can also get bug report of the whole program through statistical analysis of bug reports of all the failing runs.

Section 2 summarizes the related work; Section 3 introduces our fault localization method based on statistical differences between reduced runs; Section 4 shows comparative analysis of experimental results. Section 5 and gives the next step.

Related Work

Program runs statistical method is in essence a statistical intelligent method, how to make statistics and compare successful runs and failing runs attracted interest of many researchers [START_REF] Renieris | Fault localization with nearest neighbor queries[END_REF][START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF][START_REF] Wang | Automated path generation for software fault localization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF][START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF]. According to different type of comparison techniques, these methods can be divided into distance measures-based method and feature-based statistical approach.

Distance measures-based method is to find a successful run which is the most similar to the failing run through a certain distance measure technique, and then calculate the difference between it and failing run for fault location. A typical work of these methods is that Renieris [START_REF] Renieris | Fault localization with nearest neighbor queries[END_REF] et al proposed the "Nearest Neighbor" and Wang [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF] et al proposed calculating difference of two runs through alignment based-on control flow information.

In contrast, feature-based statistical approach locates fault-relevant statements (or faulty statements directly) by comparing the statistical information of program elements in these two kinds of run. Such program elements can be statements [START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF] or predicates [START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liblit | Bug isolation via remote program sampling[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF]. Tarantula [START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF] statistics the frequency of every statement occurs in failing runs and successful runs, and by analyzing them to get and rank the suspicious statements. Predicate-based statistical techniques, such as CBI [START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liblit | Bug isolation via remote program sampling[END_REF] and SOBER [START_REF] Liu | Sober: statistical model-based bug localization[END_REF], locate the program predicates related to faults. CBI [START_REF] Liblit | Scalable statistical bug isolation[END_REF][START_REF] Liblit | Bug isolation via remote program sampling[END_REF] measures the increase from the probability of a predicate to be evaluated to be true in all failed runs to that in all the runs, This increase is used as the ranking score, which indicates how much the predicate is related to a fault. SOBER [START_REF] Liu | Sober: statistical model-based bug localization[END_REF] defines evaluation bias to estimate the chance that a predicate is evaluated to be true in each run. In brief, CBI and SOBER use similar kinds of statistical mean comparison.

3

Fault Localization With Reduced Run

Clustering Execution Paths

Next we define execution paths and their clustering. Definition 1. (Event) An event denotes an execution of a statement of a given program. As [START_REF] Liblit | Scalable statistical bug isolation[END_REF] we use the line number i of a statement to label the event e i associated with the statement, and the statement is thus denoted as stmt(e i ). A test case that executes a given program provides the values of input variables and the expected values of output variables.

Definition 2. (Run)

A run  of a program is a sequence of events <e 0 , e 1 , e 2 ,…, e n-1 > that are sequentially associated with statements executed by a given test case. Moreover, we denote the ith event of a run π as e i π and number of events of  as ||. Consider the program in Fig. 1, we can get successful and failing runs showed in Table 1 by executing a set of test cases, where the TC inputs row sequentially put down the input values of variables x, y, and z within test cases (e.g., "2,3,4" means a test case having inputs {x=2, y=3, z=4}), while the output row describes that fact whether the TC inputs leads to a successful run or a failing run (The indexed Succ for successful runs, and Fail for failing runs), and program runs are put sequentially in columns below the corresponding test cases. For example, the TC inputs "2,3,4" yields a run <6 1 , 8 2 , 9 3 , 15 4 , 17 5 , 19 6 , 21 7 >, where the subscripts denotes the event indices (e.g., 8 2 means that the statement on line 8 is executed once and is the second event of this run). Note that some cells in program runs are empty because we align statements for highlighting the difference in program runs. Definition 3. (Dependence): Given a program P, let a variable v 1 be defined in statement S 1 and is used in statement S 2 defining value of a variable v 2 . v 2 is dependent on v 1 if changes to v 1 in the execution will is likely to influence the definition and value of v 2 . Equivalently, we also say statement S 2 is dependent on S 1 . Moreover, v 2 is control dependent on v 1 , if v 1 occurs in a condition the truth of which controls the execution from S 1 to S 2 ; otherwise, we say that variable v 2 is data dependent on v 1 . Equivalently, we raise control/data dependence to the level of statements [START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF] ). Moreover, we use dep(e i π , ) to denote the events on which e i π is dynamically control dependent in a run .

Algorithm I: delta( 1 ,  2 ) input: two program
Definition 5. (Alignment [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF]) For any pair of events e and e′ (e in run  and e′ in run ′), we say e and e′ are aligned, denoted as align(e, e′)=true, iff (1) stmt(e)=stmt(e′), and ( 2) either e and e′ are the first events appearing in  and ′ respectively, or align(dep(e, ), dep(e′, ′))=true. Definition 6. (Alignment Vector) Given two program runs π′ and π, a alignment vector, denoted as (π, π′), is an array of marks for all events in π′ and π, obtained as follows: for each event e of π, mark is 1 if there exists e′ of π′ which can aligned to e, otherwise mark is 0.

By definition 6, the following corollary is obvious:

Corollary 1. Given two program runs π′ and π, (π, π′)= (π′, π) holds。

We use Algorithm I to compute alignment vectors, where predicate alignExist(i, π 1 , π 2 ) is used to check whether there exists an event of π 2 that can be aligned to some e i of π 1 , and function alignIndex(i, π 1 , π 2 ) is used to get the index of such event of π 2 . Moreover, function .add(i) adds a mark i in the alignment vector . [START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF] where fail i (4,2,3) is abbreviated for the ith failing run with a TC inputs (4,2,3)). In the same manner, eight successful runs are divided into three classes:

(fail0,fail0)= <1,1,1,1,1,1,1>, (2) (fail0,fail1)= <1,1,1,1,1,0,1>, (3) (fail0, fail2)= <1,1,1,1,
S I :succ 0 (2,3,4), succ 4 (5,7,8); S II :succ 1 (2,4,3), succ 3 (3,4,2), succ 5 [START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF][START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF][START_REF] Liu | Sober: statistical model-based bug localization[END_REF] and succ 7 [START_REF] Liu | Sober: statistical model-based bug localization[END_REF][START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF]; S III : succ 2 (3,2,4), succ 6 (7,5,8). >, such that (1) each event e in diff(,′) is a branch event occurrence drawn from run , (2) the events in diff(,′) appear in the same order as in , that is, for all 1≤j<k,i j <i j +1 (event e ij π appears before event e ij+1 π in ), [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF] for each e in diff(,′), there exists another branch occurrence e′ in run ′ such that align(e,e′)=true (i.e. e and e′ can be aligned). Furthermore, the outcome of e in  is different from the outcome of e′ in ′. (4) All events in  satisfying (1) and ( 2) are included in diff(,′). (5) As a special case, if execution runs  and ′ have the same control flow, then we define diff(,′)=<e 0 π >. By Definition 8, the statements in the difference are all branching statements from which two runs separately go to different branches. So if there are two adjacent values in the vector such that the former is a "1" and the latter a "0", then the event that "1" represents is added to the difference. In Algorithm II, Note that .get(i) returns the ith bit of the alignment vector , while .add(i) adds a mark i in the alignment matrix . Regarding our running example, differences are summarized in Table2 and Table 3:

Difference

Fault Localization

For each predicate p in a program P, LIBLIT [START_REF] Liblit | Scalable statistical bug isolation[END_REF] estimates two conditional probabilities:(1)Pr1=Pr(P fails| p is ever observed);(2)Pr2=Pr(P fails| p is ever observed as true). LIBIT then thinks the difference Pr2Pr1 as an indicator of how relevant p is to the fault. So we conclude that difference between failing runs and successful runs indicates that faults are more likely to appear in these positions, while difference between failing runs and failing runs indicates that the probability that fault appears in these position is very small, that is, the predicate being true can not increase the possibility of failure. Let P(e,π,S) be the probability that event e of run π appears in the difference between π and successful runs, and P(e,π,F) between π and failing runs. Then for any event e in a failing run π, we define the following Score(e,π) to indicate the suspiciousness of statements stmt(e). the formulas are as follows:

( , , )

times that e appear in difference between and successful runs the total number of events in difference between and successful runs

P e S     (2)
( , , )

times that e appear in difference between and failing runs the total number of events in difference between and failing runs

P e F     (3) ( , , ) ( , , ), ( , , ) 0 ( , )
,

P e S P e F P e F Score e otherwise           (4)
Given Table 1, Table 2 and ( 2), (3), and (4), we have: P(9 3 ,F I ,S)=2/3, P(9 3 ,F I ,F)=0, and Score(9 3 ,F I )=; P(17 5 ,F I ,S)=1/3, P(17 5 ,F I ,F)=1, and Score(17 5 , FI )=1/3; P(9 3 ,F II ,S)=1/2, P(9 3 ,F II ,F)=0, and Score(9 3 ,F II )=; P(17 5 ,F II ,S)=1/2, P(17 5 ,F II ,F)=1, and Score(17 5 ,F II )=1/2. Then we get the rank score of suspiciousness that each event in the differences is the real cause of each failure. We summarize them in Table 4, where (2/3) and the like denote the case of P(e,,F)=0, P(e,,S)=2/3, in which we compare suspiciousness by P(e,,S). If the rank score of two events are equal, we consider event that appears later has a larger suspiciousness score and will be ranked in the top [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF].Assume a certain statement is responsible for a failure; it may be executed several times and appears several times in the differences. So we compute ranking scores of each statement for each failure as [START_REF] Jones | Empirical evaluation of the tarantula automatic faultlocalization technique[END_REF]. So by (5), we get Table 5 of suspicious statements for each class of failures in our running example.

( , ) ( , ) ( ( )
F I (2/3) 1/3 F II (1/2) 1/2
We can now locate faults through the analysis of Table 5. Each row of the table represents the ranking scores of each statement for a certain failing run, if there are some statements the value of which is significantly larger than others, the statements are the real cause of such failures. Similarly each column represents the ranking scores of a certain statement for each failing run,, the relevant statements with significantly larger values are the cause of these failures. As for our example, statement 9 has a large effect on both F I and FⅡ, indicating that statement 9 leads to the failure with respect to these two classes of failing runs. So we use the following equation to compute the suspiciousness of any statement.

Experiment Results

We chose four middle-sized programs with branches, by manually injecting different errors (some programs injected with an error and some with two errors). In this way, we get 43 buggy programs and 79 failing runs. Table 6 shows the description of buggy programs generated from the original four programs. To evaluate our algorithm and compared it with Wang's, we consider three cases for each failing run: (1) case "1": for ours the largest suspiciousness statement or for Wang's the whole bug report exactly indicates the actual fault position, (2) case "-1": for two methods the actual faulty location is missed in the bug report, and (3) case "0": for ours the most suspicious statement is not the actual fault position, but the actual fault location is under suspicion. Statistical results obtained from experiments are shown in Fig. 2, where the vertical axis is the number of failing runs.

Figure.2. Comparison of the bug report of two methods

Wang's algorithm returns the smallest difference as the bug report, ignoring that the position of the fault may not be in this minimal difference but in other differences. In contrast, we use statistical methods taking into account all the circumstances, and also rank the suspicious statements according to suspiciousness. The faulty statements are bound to occur in our report unless they do not appear in any of the differences.

For programs with multiple bugs, we use the evaluation criteria in [START_REF] Guo | Accurately choosing execution runs for software fault localization[END_REF] to calculate the pgm_score(P):instead of using each failing run and successful run, we use each class of failing runs and each class of successful runs. The pgm_score(P) measures the percentage of code that can be ignored for debugging, the algorithm calculates the score of the most suspicious statement and statement with its suspicious in the second place, and finally calculate pgm_score(P), as shown in Table 7.

According to table 7, when there are two errors in the program, the most suspicious statement in our bug report has higher score than that by Wang's algorithm. The second suspicious statement also has a higher score, indicating that it is also likely to be the actual fault position, that is, there are more likely two errors in the program It can be seen from Fig. 3, the time Wang spending grows exponentially as the redundant data increases since it spends much on calculating differences between each failing run and each successful run. While our algorithm eliminates redundancy by clustering all the runs, and the next difference computation is less time-consuming, so the time grows slightly. Moreover, we not only consider the suspiciousness for each failing run, but also that for the entire program through statistical analysis of all failing runs and rank statements in bug report based on their overall suspiciousness.

Conclusion

In this paper, we present an improved fault localization method using a clustering method to eliminate the path redundancy first and then by statistical analysis of differences between classes of runs to get and rank suspicious statements. Experimental results show the great improvement in terms of efficiency and accuracy in fault localization. The next step we will consider to rank all suspicious statements and do further experimental study of our techniques running against large software.

Table 1 .

 1 Test cases and execution runs of program in fig.1

	output Succ0 Succ1 Succ2 Succ3 Fail0 Fail1 Succ4 Succ5 Succ6 Succ7 Fail2 Fail3
	TCinputs 2,3,4 2,4,3 3,2,4 3,4,2 4,2,3 4,3,2 5,7,8 5,8,7 7,5,8 7,8,5 8,5,7 8,7,5
	Program runs	6 1 8 2 9 3 15 4 17 5 19 6 21 7	6 1 8 2 9 3 15 4 17 5 21 6	6 1 8 2 9 3 11 4 17 5 19 6 21 7	6 1 8 2 9 3 15 4 17 5 17 5 17 5 6 1 6 1 8 2 8 2 9 3 9 3 11 4 11 4 19 6 21 6 21 7 21 6	6 1 8 2 9 3 15 4 17 5 19 6 21 7	6 1 8 2 9 3 15 4 17 5 21 6	6 1 8 2 9 3 11 4 17 5 19 6 21 7	6 1 8 2 9 3 15 4 17 5 17 5 17 5 6 1 6 1 8 2 8 2 9 3 9 3 11 4 11 4 19 6 21 6 21 7 21 6

  runs π 1 and π 2

	output: an alignment vector (π 1 , π 2 )
	1. t=1; j=0; //temporarily store the index
		//of events of  2
	2. ={}; //the return result
	3. outer:
	4. for (i=1; i<=|π 1 |; i++){
	5.	if alignExist(i, π 1 , π 2 ) {
	6.	j=alignIndex(i, π 1 , π 2 );
	7.	k=j  temp2  1;
	8.	for (n=0; n<k; n++) .add(0);
	9.	.add(1); //add a "1"
	10.	} else {//otherwise add a "0"
	11.	.add(0); }
	12. t=j;
	13. continue outer;
	13. }
	14. return ;

  .

	Definition 4. (Dynamic Control Dependence) Given a program run π, an event π is dynamically control dependent on another event e j e i π if e j π is the nearest event coming before e i π in  such that stmt(e i π ) is control dependent on stmt(e j π

  Definition 7. (Alignment Matrix) An alignment matrix of a program run  is a matrix composed of all alignment vectors (π, π′) such that π′ is any run distinguish from π, and 0s are appended for padding if the length of alignment vectors are not constant.By TableIand definition 6, we have[START_REF] Jones | Visualization of test information to assist fault localization[END_REF]

  1,1,1>, and (4)

	(fail0,fail3)= <1,1,1,1,1,0,1>. With such an alignment matrix, we cluster the matrix
	by column and thus four failing runs in Table I can be divided into two classes:
	F I : fail 0 (4,2,3), fail 2 (8,5,7); F II : fail 1 (4,3,2), fail 3
	Algorithm II: diff( 1 ,  2 )
	input: two program runs π 1 and π 2
	output: a set of event indices
	representing
		difference between π 1 and π 2
	1. ={};
	2. =delta( 1 ,  2 );
	3. n=0; //to store the total number
		//of 0s coming before 1
	4. for (i=0; i <||  1; i++) {
	5.	if (.get(i)=0) n++;
	6.	if (.get(i)=1.get(i+1)=0) {
	//if the ith event can be aligned, while
	//the next event can't be aligned
	7.	.add(i1n||); }
	8. }

Table 2 .

 2 Difference between each class of failing runs and successful runs

		S I	S II	S III
	F I	9 3	9 3 ,17 5	null
	F II	9 3 ,17 5	9 3	17 5

Table 3 .

 3 Difference between each class of failing runs

		F I	F II
	F I	null	17 5
	F II	17 5	null

Table 4 .

 4 Score of suspicious events for failures

	Score of suspicious events of F I
	F I	9 3 (2/3)	17 5 1/3

Score of suspicious events of F II

	F I F II	9 3 (1/2)	17 5 1/2

Table . 5

 . . Score of suspicious statements for failures

	s	Statement	9	17

Table 6 .

 6 Description of experiment data

	Progr ams P1	nested level branches 0	of	having one bug 5 a	Buggy Versions Kinds of having failing runs two bugs 14 3	Kinds of failing runs 11
	P2	1		4	4	4	7
	P3	2		8	9	8	16
	P4	3		5	5	7	13

Table 7 .

 7 Comparison of two methods when there are two errors in the program

					Figure.3 Comparison of run time
					of two methods
	Score 0.8-1.0 0.7-0.79 0.6-0.69 0.5-0.59	Wang algorithm 13.3 60.0 6.7 13.3	The first place of our algorithm 33.3 26.7 6.7 26.7	The second place of our algorithm 46.7 6.7 20.0 6.7	time(ms)	1000 900 800 0 100 200 300 400 500 600 700	1	2	3	4	repeat counts 5 6	7	8	9	1 0
	0-0.49	6.7	6.7	20.0					Wang Algorithm	Our Algorithm

The corresponding Author. This work was supported by National Natural Science Foundation of China (60775028), Dalian Science &Technology Program (2007A14GX042), IT Industry Development of Jilin Province.


Dmitry Devetyarov
email: dmitry@cs.rhul.ac.uk

Ilia Nouretdinov

Prediction with Confidence Based on a Random Forest Classifier

Keywords: Conformal predictor, confidence machine, region prediction, random forest

Conformal predictors represent a new flexible framework that outputs region predictions with a guaranteed error rate. Efficiency of such predictions depends on the nonconformity measure that underlies the predictor. In this work we designed new nonconformity measures based on a random forest classifier. Experiments demonstrate that proposed conformal predictors are more efficient than current benchmarks on noisy mass spectrometry data (and at least as efficient on other type of data) while maintaining the property of validity: they output fewer multiple predictions, and the ratio of mistakes does not exceed the preset level. When forced to produce singleton predictions, the designed conformal predictors are at least as accurate as the benchmarks and sometimes significantly outperform them.

Introduction

The new framework of conformal prediction introduced in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] allows us to output region predictions (a set of predicted labels) with the guaranteed error rate under a simple statistical assumption (this property is called validity), as opposed to point predictions when we always produce singleton predictions but the error rate is not guaranteed.

Having a guaranteed error rate, we may still obtain multiple region predictions, that is, predictions that comprise more than one label. Multiple predictions are not mistakes: they indicate that there was no sufficient information provided for predicting one label. The question is how big these region predictions are, and our aim is to decrease a number of multiple predictions. The ability of conformal predictors to produce predictions as certain as possible is called efficiency.

Most of known machine learning algorithms can be used as an underlying algorithm in a conformal predictor. Efficiency of conformal predictors is usually in line with the accuracy of the underlying algorithm and therefore varies across the range of underlying algorithms and also depends on the type of data analysed. For this reason, we are looking for new nonconformity measures that could result in efficient predictions. So far various nonconformity measures have been designed. In this paper we propose new conformal predictors based on a random forest. We expect the predictors to inherit random forest advantages and to maintain the property of validity.

Although the main aim of this paper is to elaborate and analyse new methodology of providing region predictions, we can force conformal predictors to output one prediction instead of multiple predictions. Such approach allowed us to compare designed conformal predictors with machine learning methods (random forest, in particular). But one should bare in mind that in this case we do not have the advantages of conformal predictors as region predictors: there is no guaranteed validity.

Outline of Conformal Prediction

The framework of conformal prediction is described in detail in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. Conformal predictors are based on the only assumption about the data generating mechanism: all the examples have been generated independently by some probability distribution (the i.i.d. assumption).

Let us assume that we are given a training set of examples (x 1 , y 1 ), . . . , (x n-1 , y n-1 ), where x i ∈ X is a vector of attributes and y i ∈ Y is a label out of a finite set of possible labels (classes), and our goal is to predict the classification y n for remaining example x n . A combination of an example and a label

z i = (x i , y i ) ∈ Z = X × Y is called an object.
A nonconformity measure is a set of measurable mappings 1) is the set of all bags (multisets) of elements of Z of size n -1. For each possible label y, we consider the hypothesis y n = y and the nonconformity measure assigns some values α i (nonconformity scores) to every example in the sequence {z i , i = 1, . . . , n} including a new example and evaluates 'nonconformity' α i := A n ( z 1 , . . . , z i-1 , z i+1 , . . . , z n , z i ), i = 1, . . . , n between a set and its element ( . . . denotes a multiset).

{A n : n ∈ N } of the type A n : Z (n-1) × Z → (-∞, +∞], where Z (n-
For each hypothesis y n = y, we compare α n to the other α i s and calculate p(y) = |{i = 1, . . . , n : α i ≥ α n }|/n -the p-value associated with the possible label y for x n . Thus, we can compliment each label with a p-value that shows how well a new example with this label conforms with the rest of the sequence.

The conformal predictor determined by the nonconformity measure A n , n ∈ N and a significance level is then defined as a measurable function Γ :

Z * × X × (0, 1) → 2 Y (2 Y is a set of all subsets of Y ) such that the prediction set Γ ( ) (x 1 , y 1 , . . . , x n-1 , y n-1 , x n ) is
defined as the set of all labels y ∈ Y such that p n > . Thus, for any finite sequence of examples with labels, (x 1 , y 1 , . . . , x n-1 , y n-1 ), a new unlabelled object x n and a significance level , the conformal predictor outputs a region prediction Γ ( ) -a set of possible labels for a new object.

Conformal predictors defined above are valid : in the long run the frequency of errors made by a conformal predictor (that is, cases when prediction set Γ does not contain a true label) does not exceed subject to the i.i.d. assumption. Strictly speaking, for any exchangeable probability distribution P on Z ∞ (Z ∞ is the set of all infinite sequences of elements of Z) and any significance level , lim sup

n→∞ n i=1 err n (Γ ) n ≤ (1)
with probability one, where err n (Γ ) is equal to 1 when the prediction set Γ does not contain a real label y n , and 0 otherwise. The property of validity is theoretically proven in the on-line mode and empirically confirmed in the offline mode [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF].

Random Forests

In this work we consider the type of random forests described in [START_REF] Breiman | Random Forests[END_REF]. Theoretical results [START_REF] Breiman | Random Forests[END_REF] demonstrate that random forests do not overfit when more trees are added. They also empirically proved to have the following advantages ([4], [START_REF] Breiman | Random Forests[END_REF]): random forests produce high accuracy for many data sets; they can process data with a large number of features where each feature is weak, that is, carries a small amount of information; they are relatively robust to mixed variable types, missing data, outliers and noisy data; constructing random forests is relatively fast (faster than bagging and boosting).

In brief, a random forest is a classifier that consists of decision trees, each of which provides a vote for a certain class. Combining a large number of trees in a random forest can lead to more reliable predictions, while single decision tree may overfit the data.

Nonconformity Measures Based on Random Forests

In this paper we designed nonconformity measures based on random forests.

We will use the following notation: suppose we are given a bag (x 1 , y 1 ), (x 2 , y 2 ), . . . ,(x m , y m ) , (x i , y i ) ∈ Z, and we need to define a nonconformity measure A(x, y) = A( (x 1 , y 1 ), (x 2 , y 2 ), . . . , (x m , y m ) ; (x, y)). Alternatively, we can define a conformity measure B(x, y) = 1 -A(x, y) when it is more intuitive.

The nonconformity or conformity measures we propose are the following: 1. A random forest is constructed from a training set {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x m , y m )}. The conformity score of a new example (x, y) is then equal to the percentage of correct predictions given for x by decision trees.

2. Conformity measure 1 is the most natural one, however, it is computationally inefficient: when considering example N + 1 we have to construct (N + 1)L random forests, where L is the number of labels. We will therefore use another conformity measure, which will require only one random forest when making a prediction for a new object. The random forest is grown for the union of a bag (x 1 , y 1 ), (x 2 , y 2 ), . . . , (x m , y m ) and a new example (x, y). Since for each decision tree, the training set is a bootstrap sample, a new example is not included in this training set in about one third of decision trees. For each (x, y) we aggregate the votes for this example only of those decision trees where this example is out-of-bag (not in the training set for the tree). The conformity score is then equal to the proportion of correct votes for (x, y) among these trees.

3. This nonconformity measure was proposed by Huazhen Wang and Fan Yang in our personal communication. It is based on random forest proximities P (i, j), i, j = 1, . . . , m + 1, which provide a measure of how close to each other two objects are regardless of their labels and are calculated as a ratio of trees, running through which objects i and j land at the same terminal node. We construct a random forest for the union of a bag (x 1 , y 1 ), (x 2 , y 2 ), . . . , (x m , y m ) and a new example (x, y) and form the corresponding (m+1)×(m+1) matrix of proximities for objects

i s and j s are the numbers of examples with s-st greatest value of proximity with example (x, y) among examples labelled with the same label y and among all the other examples, respectively.

Experiments

The designed nonconformity measures were implemented and applied to different data sets.

Data

In our experiments we used six proteomic data sets, two medical non-proteomic data sets and two non-medical data sets. Proteomic data sets comprise: ovarian cancer data from the UKOPS trial [START_REF] Timms | Peptides Generated Ex Vivo from Abundant Serum Proteins by Tumour-Specific Txopeptidases are Not Useful Biomarkers in Ovarian Cancer[END_REF]; ovarian cancer (OC ), breast cancer (BC ) and heart disease (HD) samples collected in the UKCTOCS trial1 ; Competition data provided by the Leiden Clinical Mass Spectrometry Proteomic Diagnosis Competition and preprocessed as described in [START_REF] Gammerman | Clinical Mass Spectrometry Proteomic Diagnosis by Conformal Predictors[END_REF]; Ciphergen's 7 biomarkers of UKOPS data [START_REF] Nouretdinov | Data Analysis of 7 Biomarkers[END_REF]. The other medical data sets are shortened Abdominal pain data [START_REF] Gammerman | Bayesian Diagnostic Probabilities without Assuming Independence of Symptoms[END_REF], which comprises 33 symptoms of acute abdominal pain, and Microarray data of lung cancer, colon cancer and breast cancer patients provided in the ICMLA 2009 Challenge [START_REF] Nouretdinov | Application of Inductive Confidence Machine to ICMLA Competition Data[END_REF]. Sonar and Iris non-medical data sets were taken from the University of California, Irvine (UCI) Machine Learning Repository.

Implemented Conformal Predictors

We implemented non-conformity measures 2 (referred to as CP-RF ) and 3 (referred to as CP-RF-kNN, where k is the number of nearest neighbours). These conformal predictors were compared with benchmark predictors based on the kNN algorithm [3] (we will denote them CP-kNN, where k is the number of nearest neighbours) and SVM with different kernels [START_REF] Gammerman | Learning by Transduction[END_REF] (denoted as CP-SVM (kernel, parameter )).

The experiments were carried out in two settings: off-line in leave-one-out (LOO) procedure, to show the usage of conformal predictors as conventional classifiers, and on-line, to demonstrate the advantages of region predictions.

We used the following parameters for random forest construction: the number of trees 1000, the number of features selected at each node to split equals a square root of the number of features. These values are recommended in [START_REF] Breiman | Random Forests[END_REF], where it is theoretically proven that the results converge when we increase the number of trees in random forests and it is empirically shown that the results are insensitive to the number of features selected at each node.

Results

For each significance level > 0, conformal predictors output a set of labels with p-values greater than . Thus, the predictor may output no label (we call this an empty prediction), one label (certain prediction) or more than one label (multiple prediction).

Firstly, the designed conformal predictors proved to be valid, that is, for a given significance level > 0 the rate of erroneous predictions, that is, predictions not containing an actual label, is close to .

The example of the erroneous prediction dynamics is shown in Figure 1a. The figure demonstrates validity of the CP-RF-1NN applied to the Microarray data for significance levels = 5% and 10% : solid lines, which represent the actual number of errors, are close to dotted lines, which demonstrate the expected number of errors for different significance levels.

Figure 1b demonstrates the dynamics of efficiency characteristics at significance level of 10% of the CP-RF-1NN applied to the Microarray data in the on-line mode. The characteristics shown are the number of multiple predictions, the number of certain predictions and the number of empty predictions. The figure demonstrates that while the number of analysed examples is low, they do not carry enough information to make certain predictions without losing validity. But starting from example 50, we have accumulated enough information so that multiple predictions cease to occur and most of prediction regions contain exactly one label, which is in most cases correct. The dynamics on the plot also conforms with the empirical fact established in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] that when multiple predictions disappear, empty predictions start to occur.

As mentioned before, all implemented conformal predictors have a theoretically proven property of validity, and the general aim is to design a nonconformity measure that could improve efficiency, that is, make the algorithm output as few multiple predictions and as many empty predictions as possible. Comparison of efficiency of CP-RF and CP-RF-kNN with the benchmarks demonstrated that CP-RF and CP-RF-kNN produce at least as few multiple predictions and as many correct certain and empty predictions as the benchmarks, and they perform much better in terms of efficiency on all mass spectrometry data sets. This allows us to speculate that conformal predictors based on random forests benefit from the advantages of the underlying algorithm and perform well on noisy data and data with a lot of weak inputs. Table 1 summarizes the multiple prediction rate for different conformal predictors. Conformal predictors have been developed to provide region predictions with the preset error rate. However, in order to compare them with bare predictions output by conventional machine learning methods, we can ignore the nature of conformal predictors and force them to always make a certain prediction. After assigning a p-value for each label to every object, we can predict the label with the highest p-value. This is called forced point prediction. If several labels have the highest p-value (we call this situation a tie), we make a random prediction.

Experiments demonstrated that when forced to make point predictions, conformal predictors perform similarly to random forest algorithm (see Table 2). This can be explained by the fact that each random forest is a combination of a large number of trees constructed randomly and each sample is not included in the training set for about one third of all trees in a random forest.

This implies that we can add the framework of conformal prediction to the random forest algorithm without losing in accuracy, while benefiting from conformal predictions: we can produce valid region predictions and compliment each prediction with confidence.

The results of comparison of forced point prediction accuracy of different conformal predictors (Table 2) were in line with efficiency comparison: CP-RF and CP-RF-kNN significantly outperformed other predictors on certain mass spectrometry datasets and were at least as good as the benchmarks on all data sets.

Conclusion

In this paper we worked on further development of conformal predictors, which produce region predictions that make a preset number of errors in the long run. Designed nonconformity measures based on random forests proved to be valid and efficient. First, the ratio of mistakes does not exceed the preset level. Second, CP-RF and CP-RF-kNN produce more efficient predictions on all mass spectrometry data sets and are not beaten on the other data. When forced to produce singleton predictions, conformal predictors based on random forest result in accuracy similar to random forest accuracy. The accuracy of forced point predictions output by CP-RF and CP-RF-kNN is at least as high as accuracy produced by known conformal predictors and sometimes is significantly higher. This implies that although conformal predictors are designed for producing valid region prediction, they can also be a useful tool when making singleton predictions.

  x 1 , x 2 , . . . , x m , x m+1 = x. The nonconformity measure is the ratio of the average proximity of the example with examples of other classes to the average proximity of the example to examples of the same class. In both averages we consider only proximities of those k examples that have the greatest values of proximities among examples of the same class y and among all the other examples. Strictly speaking, A(x, y) = A(x, y) -/A(x, y) + , where A(x, y) + = k s=1 P (i s , m + 1), A(x, y) -= k s=1 P (j s , m + 1) ,

Fig. 1 :

 1 Fig. 1: Validity and efficiency of CP-RF-1NN applied to the Microarray data in the on-line mode.

Table 1 :

 1 The rate of multiple predictions for significance level = 10% in the LOO mode.

	Data	CP-RF CP-RF CP-RF CP-1NN CP-5NN CP-SVM CP-SVM
		-1NN -5NN			(rbf, 5) (poly, 5)
	UKOPS	46.1% 47.0% 45.8% 74.8% 72.0%	59.2%	69.8%
	UKCTOCS OC 16.0% 16.0% 13.8% 44.6% 30.8%	38.5%	79.8%
	UKCTOCS BC 77.8% 78.4% 77.8% 80.9% 80.9%	81.5%	82.7%
	UKCTOCS HD 56.0% 58.1% 57.2% 64.7% 59.5%	66.1%	64.0%
	Competition 11.1% 18.3% 17.0% 26.8% 19.6%	32.7%	30.7%
	7 biomarkers 51.1% 55.4% 53.8% 67.0% 61.2%	97.9%	96.9%
	Abdominal pain 0.3% 1.0% 0.0%	3.0%	0.0%	0.0%	1.0%
	Microarray	0.0% 1.5% 0.3% 13.5%	3.5%	8.5%	40.4%
	Sonar	14.9% 11.1% 13.0% 13.9% 16.4%	32.2%	30.8%
	Iris	0.0% 0.0% 0.0%	0.0%	0.0%	86.7%	8.0%

Table 2 :

 2 Accuracy of forced point predictions in the LOO mode.

	Data	RF CP-RF CP-RF CP-RF CP-1NN CP-5NN CP-SVM CP-SVM
		-1NN -5NN	(rbf, 5) (poly, 5)
	UKOPS	72.6% 72.3% 71.5% 72.6% 55.1% 61.7%	66.7%	55.5%
	UKCTOCS OC 84.9% 84.8% 83.8% 84.6% 72.3% 80.6%	78.9%	77.9%
	UKCTOCS BC 66.0% 66.7% 59.0% 62.4% 50.3% 62.4%	56.2%	54.3%
	UKCTOCS HD 71.8% 72.3% 69.2% 71.4% 63.2% 67.9%	62.4%	62.1%
	Competition 83.7% 85.3% 83.7% 83.3% 82.0% 84.6%	86.3%	87.6%
	7 biomarkers 74.6% 74.8% 72.2% 73.9% 64.5% 73.7%	60.9%	59.0%
	Abdominal pain 91.7% 92.7% 91.8% 91.5% 88.0% 92.2%	91.7%	90.7%
	Microarray 92.0% 91.3% 92.8% 91.4% 86.1% 89.4%	88.3%	89.5%
	Sonar	85.1% 84.6% 88.7% 85.6% 86.3% 82.9%	84.6%	85.3%
	Iris	95.3% 94.7% 95.0% 95.3% 93.3% 97.0%	89.3%	89.7%

for more information see www.ukctocs.org.uk

Acknowledgments. We would like to thank our collaborators in MRC project ("Proteomic analysis of the human serum proteome") Ian Jacobs, Usha Menon, Rainer Cramer, John F. Timms, Ali Tiss, Jeremy Ford, Stephane Camuzeaux, Aleksandra Gentry-Maharaj for collecting the original UKCTOCS and UKOPS data and carrying out MS experiments. This work is supported by EPSRC grant EP/F002998/1 (Practical competitive prediction); MRC grant G0301107 (Proteomic analysis of the human serum proteome); VLA of DEFRA grant on Development and Application of Machine Learning Algorithms for the Analysis of Complex Veterinary Data Sets; EU FP7 grant O-PTM-Biomarkers HEALTH-2007-2.4.1-2; ASPIDA grant "Development of new methods of conformal prediction with application to medical diagnosis" from the Cyprus Research Promotion Foundation; European Union EU FP7 programme "Living with uninvited guests: comparing plant and animal responses to endocytic invasions" (ERASysBio); and MRC Award G0802594.


François Rioult

Sami Mecheri

Bruno Mantel

François Kauffmann

Nicolas Benguigui

What Can Hawk-Eye Data Reveal about Serve Performance in Tennis?

In the present study, we aim at showing how some characteristics of the serve summed up in the resulting ball trajectory can determine the efficiency of tennis serves. To that purpose, we analyzed a big set of data collected between 2003 and 2008 at international ATP, WTA and Grand Slam tournaments and corresponding to 84 tournaments, 1729 matches, 262,596 points. Using time-dependent three-dimensional ball trajectory data recorded by the automated ball tracking Hawk-Eye system, we show the relationships that exists between the characteristics of the serve kinematics and impacts on the ground on the gain of the points.

Introduction

In recent years, the development of technologies and automatic tracking systems has enabled the capture of ball trajectories during tennis matches. Since 2003, Hawk-Eye vision-based systems have provided ball tracking to assist the players when they think an error of judgments has been made by referees. This system uses a motion capture system with 10 cameras around the court and sophisticated algorithms will calculate the trajectories and impact on the ground of the tennis ball with an accuracy estimated of 3.6 mm at impact.

Despite the high-level accuracy of such tracking systems and the huge amount of kinematic data generated, the use of these systems for quantitative analysis of player performance and scientific analysis is rare and has never been performed on a very big sheer volume of data. To our knowledge, only three studies used Hawk-Eye data to analyze performance. These studies were interested in prediction of shot locations ( [START_REF] Wei | Predicting shot locations in tennis using spatiotemporal data[END_REF], data volume = matches from the Australian Open men's draw or around 10 000 points), in laterality effect on ball distribution ( [START_REF] Loffing | The serve in professional men's tennis: Effects of players' handedness[END_REF], data volume = 32 matches or 4744 points) and in on-court position effect on groundstroke anticipation ( [START_REF] Loffing | On-court position influences skilled tennis players' anticipation of shot outcome[END_REF], data volume = 38 matches, number of points unspecified).

In the present study, we aim at showing how some characteristics of the serve summed up in the resulting ball trajectory can determine the efficiency of tennis serves. To that purpose, we analyzed a big set of data collected between 2003 and 2008 international at ATP, WTA and Grand Slam tournaments and corresponding to 84 tournaments, 1729 matches, 262,596 points.

The influence of factors such as serve speed, serve location, court-surface and men/women differences on the winning-point rate was assessed in order to provide an extensive insight into efficient serve tendencies in world-class tennis. The positions of serves' impact were also examined in order to provide an accurate description of the serves performed by world-class players during matches. Since the present work is the first to exploit large-scale Hawk-Eye data, a subsidiary objective in these analyses was to demonstrate our method as reliable to analyze serving match strategies by confronting our findings to knowledge emanating from tennis performance analysis studies [START_REF] Gillet | A notational analysis of elite tennis serve and serve-return strategies on slow surface[END_REF][START_REF] O'donoghue | A notational analysis of elite tennis strategy[END_REF].

We also focused on the unexplored question of the magnus effect intensity in serve trajectories. The spinning of the tennis ball was characterized in the present study directly from kinematic data by the ball axis of rotation and the speed of rotation around this axis. Specifically, the lift coefficient (as an indicator of spin intensity) and the ball axis of rotation (as an indicator of spin nature) were analyzed.

Data description

The data analyzed in the present research were made available by the company Hawk-Eye Innovations in the context of a publicly funded research project (Ten-nisServer, ANR-06-BLAN-0413) in which one of us was involved in 2006-2009. For the moment being, the data are not publicly available.

The final stages of the most famous tournaments of the ATP and WTA circuits between 2003 and 2008 are covered by the data. 40 Hz trajectory of the ball and XML information about the points are available. Each file is named after the number of the set, the number of the game, the index of the point, the index of the serve (first or second, there is no file in case of double fault), and the time of the point.

For each point, the XML file gathers the following information (see Figure 1 for an excerpt):

1. the header gives overall information about the point: the server, the receiver, the player who is located on the positive part of the court, the class of the serve (0 for an ace, 1 for a classical one, 2 for a winning serve), the scorer of the point (1 if he/she is the server, -1 otherwise), the duration of the point (in seconds), and the score in the game at the start of the point; 2. the precise coordinates of the serve: who serves, the initial speed, the final speed, the coordinates of the initial impact (at t = 0), the coordinates of the bounce; 3. the precise coordinates of each shot.

After cleaning the data, there remains 75,587 points for the women and 187,009 for the men (total: 262,596 points).

<?xml version="1.0" encoding="UTF-8" standalone="no" ?> <point valid="true"> <hawkeye_header> <xmldate d="Data"/> <server p="CLIJSTERS"/> <receiver p="HENINHARDENNE"/> <positive p="HENINHARDENNE"/> <serve_class c="1"/> <scorer s="-1"/> <PointDuration w="6.26786"/> <score_raw s="1 0"/> </hawkeye_header> <serve name="CLIJSTERS" player="1" speed="46.46" speedEnd="31.87"> <coord t="0" x="1.47" y="-11.89" z="2.70"/> <coord bounce="true" t="0.49025" x="-3.05" y="6.17" z="0.033"/> </serve> <shot speed="31.3742" speedEnd="19.8407"> <coord t="0.8" x="-4.15587" y="11.766" z="1.06128"/> <coord bounce="true" t="1.57191" x="0.24" y="-6.26" z="0.033"/> </shot> ... </point> Fig. 1: Sample of XML data.

Linear magnus model

In this section, we aimed to model the kinetics of a spinning tennis ball by estimating unknown parameters from reconstructed trajectories, using the R software [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Our analysis revealed that the Hawk-Eye reconstructed trajectories are using a third degree polynomial in each components (x, y, z).

In contrast with previous studies which obtained spin rates by manually counting the number of revolutions of the ball from high-speed video cameras recordings (e.g., maximum serve spin rates values of 3529 rpm in Wimbledon qualifications reported by [START_REF] Kelley | Ball spin generation at the 2007 Wimbledon qualifying tournament (P110)[END_REF] and of 4300 rpm in Davis Cup reported by [START_REF] Goodwill | Measuring ball spin off a tennis racket[END_REF]), we used reconstructed ball trajectories and characterized the spinning of the tennis ball by its axis of rotation ω and the speed of rotation around this axis ω.

We use the model proposed in [START_REF] Klvana | Trajectory of a Spinning Tennis Ball[END_REF] to simulate the tennis ball trajectories. The tennis ball is considered as a mass point at position X(t) = (x(t), y(t), z(t)) with mass m, diameter d and is influenced by three forces :

the weight force G = mg with g = (0, 0, -g)

-the drag force D = -D L (v, ω) V v with D L (v, ω) = C D (v, ω) 1 2 πd 2 4 ρv 2 -the magnus force M = M L (v, w) ω ω ∧ V v with M L (v, ω) = C L (v, ω) 1 2 πd 2 4 ρv 2
We introduce physical characteristics of a reference tennis ball and atmospheric conditions:

a reference diameter d 0 = 0.067m a reference density of the air ρ 0 = 1.29 a reference mass m 0 = 0.0577kg If we write α = πd 2 ρ 8m and α 0 = πd 2 0 ρ0 8m0 then the equation becomes:

d 2 X(t) dt 2 + 9.81g = C G g -C D α α 0 (α 0 vV ) + C L α α 0 ω ω ∧ (α 0 vV ) .
The scalar coefficient C G may be interpreted as a variation between normal gravity 9.81m/s 2 and the observed gravity during the match. Variability of this coefficient may be due to latitude variation ±0.03m/s 2 , gravitational anomalies, sampling frequency or model errors.

If we define the modified drag coefficient as

C d = C D α α0
and the modified lift coefficient as C L = C L α α0 then the equation becomes:

d 2 X(t) dt 2 + 9.81g = C G g + C D (-α 0 vV ) + C L ω ω ∧ (α 0 vV ) (1)
The main assumptions of this magnus linear model is that the modified drag and lift coefficients are constant throughout an arc.

In Equation 1, the vectors d 2 X(t) Alam [START_REF] Alam | The aerodynamic forces on a series of tennis balls[END_REF] has estimated the drag coefficient C D in the absence of any spin to lie between 0.5 and 1.2 for the tennis ball. At lower velocity, the mean value was 0.90 ± 0.15, whereas at higher velocity the mean value was 0.6 ± 0.025.

Goodwill [START_REF] Goodwill | Aerodynamics of spinning and nonspinning tennis balls[END_REF] has studied the lift coefficient C L of a tennis ball in a wind tunnel in different conditions, as a function of S = d/2ω v . Drag coefficients were varying from 0.65 ± 0.01 for low value of S i.e. S ≤ 0.3 and raised to 0.69 ± 0.01 for higher S values. [START_REF] Goodwill | Aerodynamics of spinning and nonspinning tennis balls[END_REF][START_REF] Brody | The Physics and Technology of Tennis[END_REF][START_REF] Cross | Measurements of drag and lift on tennis balls in flight[END_REF] found lift coefficient from 0.02 to 0.3.

We found that 80 % of points' trajectories had a global R 2 greater than 0.97, meaning that for this subset of points, the linear combination of these three estimated components

C G g, C D (α 0 vV ), C L ω ω ∧ (α 0 vV ) provided a good ap- proximation of d 2 X(t) dt 2 + 9.81g.

Results

Winning probability for server

The results of Table 1 shed light on the fact that the serve is a redoubtable shot for winning points in tennis. It provided servers with the opportunity to accumulate a high percentage of winning points, particularly from the first serve (69.46%4 ). This advantage of the server over the receiver confirms the results of [START_REF] O'donoghue | A notational analysis of elite tennis strategy[END_REF][START_REF] Unierzyski | Comparison of tactical solutions and game patterns in the finals of two grand slam tournaments in tennis[END_REF] as well as [START_REF] Gillet | A notational analysis of elite tennis serve and serve-return strategies on slow surface[END_REF] who reported 67.3% wins and 53.8% wins on second serves on clay (66.28% and 52.24% in the present study). Unsurprisingly, the court surface also had a significant influence on winning rate.

Impact of the gender

The analysis (see Table 2) revealed that men won significantly more points when serving than women both on first and second serves. Other research efforts have also noted gender differences in winning percentages on serve [START_REF] O'donoghue | The impact of speed of service in Grand Slam singles tennis[END_REF][START_REF] O'donoghue | A notational analysis of elite tennis strategy[END_REF][START_REF] Verlinden | Effect of Gender and Tennis Court Surface Properties Upon Strategy in Elite Singles[END_REF]. This result could be mainly explained by the difference in speed of serves across men and women.

Impact of serve speed

The results of Figure 1 indicate a significant effect of serve speed on winning percentage on serve. These results are in agreement with the findings of [START_REF] O'donoghue | The impact of speed of service in Grand Slam singles tennis[END_REF] who have noted a significant relationship between the serve speed and the probability of winning the point. They found that, serve speed was negatively correlated with the proportion of serves that fell inside the serve box. Also, the proportion of points won when the serve was in was positively correlated with the serve speed for both the first and second serves in Grand Slam tournaments [START_REF] O'donoghue | The impact of speed of service in Grand Slam singles tennis[END_REF]. Therefore, hitting a "hard" first serve is a winning serve strategy to win a high percentage of points [START_REF] Brody | The Physics and Technology of Tennis[END_REF]. This strategy increases the time constraints on receivers by reducing the time available for executing their shot.

Impact of the number of strokes

Results of Figure 2 showed a clear negative relationship between the number of recorded shots per rally and first serve's winning percentage in both men and women. These data can be summarized in the following way: the lower the number of strokes per point, the greater the impact of serve on winning the point. This makes sense: the severe spatio-temporal constraints imposed on the receiver when facing a first serve make it difficult to restore the balance in one or two groundstrokes. Consequently, points concluded within a very low number of shots result to outcomes that tend to favor the servers. However, for second serves, the relationship between the number of recorded shots per rally and the serve winning percentage is much different since the [0,2] category is linked to a very low winning rate (44% in men and 34% in women) while other categories all share similar winning rates (around 43% in men and 45% in women). Since the spatio-temporal constraints on the receiver are considerably lower for second serves, the receiver has the opportunity to initiate a baseline rally in almost all cases. In this configuration, if the point is gained quickly, it cannot be attributed to serve impact (whose speed is very moderate) but to the quality of the return. These results nevertheless remain surprising for their magnitude, and another explanation could be problems in the sub-sets of data used to perform these analyses.

Serves' impacts positions

Fig. 4: Distribution of serves' impacts (y-axis) in service boxes for the men (left) and the women (right). Deuce court is located on the left of the server while the advantage court is on his/her right.

The distribution of serves' impacts positions along the y-axis5 (Figure 4) revealed that most first serves were directed toward the T (middle) and W (edge) locations in a very similar fashion on both deuce and advantage courts. A similar trend has already been reported in a previous work describing the serve locations of male professionals on hard courts [START_REF] Unierzyski | Comparison of tactical solutions and game patterns in the finals of two grand slam tournaments in tennis[END_REF]. The first serve is typically at pace [START_REF] Magnus | On the advantage of serving first in a tennis set: four years at Wimbledon[END_REF], and one would expect first serves to be directed more often toward the W and T locations because serving to these locations takes the ball away from the receiver, making it difficult for him/her to return (which is the main goal of first serves [START_REF] Crespo | Advanced Coaches Manual[END_REF]).

Interestingly, the distribution of serves' impacts positions along the y-axis revealed differences between deuce (mainly T) and advantage courts (mainly W) for second serves. This finding is in line with [START_REF] Gillet | A notational analysis of elite tennis serve and serve-return strategies on slow surface[END_REF] which noted that close to 95% of second serves were directed either toward the T in the deuce court (48.0%) or to the wide zone in the advantage court (46%). In other words, on second serve on both sides of the court, professional tennis players serve to the corners of the service box with a specific focus on their opponents' backhand (most of which are right-handed opponents), which is usually considered the weaker side. Our data strongly confirmed that two strategies are employed on second serves, depending on the service box being played. When serving on the deuce side, servers attempt to push back the receivers with a topspin (as demonstrated by C L values) toward the T so as to keep them behind the baseline. When serving on the advantage side, players attempt to find more angles by serving wide and with topspin (as demonstrated by C L values) to open up the court. In both cases, the server's intention is to dominate the rally from its start by exerting a territorial influence.

Lift coefficient and rotation axis

The analysis (see Figure 5) revealed clear differences in C L as a function of serve ball with C L values significantly higher for second serves. This result confirms and extends current knowledge about tennis. Indeed, [START_REF] Gillet | A notational analysis of elite tennis serve and serve-return strategies on slow surface[END_REF] reported that for first serves, the flat option was the most used (55.7%) while for second serves spin variations were massively used (99.0%). First serve spins are employed to introduce tactical variations but with parsimony since it reduces serve's speed. However, during second serve, the players' goal is to limit aggressive and offensive returns. For this reason, as reported by [START_REF] Gillet | A notational analysis of elite tennis serve and serve-return strategies on slow surface[END_REF], topspin strategy is classically used on the second serve (91.6%) to generate a shoulder or head-high and deep bounce, which prevents the receiver from executing an offensive stroke.

Conclusion

The present study has confirmed and extended knowledge about tennis duel by manipulating various performance indicators of the first stroke of each point and assessing its influence on winning-point probabilities. Having demonstrated the validity of our trajectory reconstruction method for tennis performance analysis by replicating the findings of earlier tennis performance studies ([4, 14, 15]), this method could be used to provide coaches and researchers with objective and massive information on serving performance. On this plan, the high proportion of first serves oriented to the two corners of the box revealed that players do not maximize the possibility of varying the direction of the serve. The above results could be made even more meaningful by incorporating serve variability indicators such as entropy to determine the location succession effect on serve winning rate. Future research on this plan is encouraged to disentangle the complexity of situational probability information that is integrated into decisions of expert players in serve-return.

Also, the direct method of spinning determination used in the present paper is highly valuable since it is applicable with no supplementary time costs to all players competing or just training on ball tracking equipped-courts. Obtaining spinning data by this way is desirable since it allows players and coaches to obtain accurate information about their stroke quality from matches and practice that are readily available by avoiding manual analysis which is time consuming. This approach based on 3d-ball tracking data not only offers further works for serve or serve-return performance, but might also help to add further knowledge about players' fitness level and ground-stroke quality.

dt 2 + 9 .

 29 81g and α 0 vV can be estimated with the model for speed and acceleration. The four unknown coefficients are C G , C D , C L ωx ω , C L ωy ω , C L ωz ω , appear linearly in the equation and therefore may be estimated with a linear model. Modified drag and lift coefficients C D , C L depend on properties of the roughness of the ball's surface, on velocity and on spinning. For a tennis ball which has the characteristic of the reference ball we have C D = C D and C L = C L . The factor α α0 is a correction factor which only depends on the cross sectional area and the mass of the tennis ball in comparison to a reference tennis ball.

Fig. 2 :

 2 Fig. 2: Probability to win a point w.r.t. the speed of the serve (lose on top, win on bottom).

Fig. 3 :

 3 Fig. 3: Probability to win according to the number of strokes.

Fig. 5 :

 5 Fig. 5: Distribution of φ as a function of serve and modified lift coefficient C L for right handedness (φ is the angle of the rotation vector in the yz-plane, φ = 180 o means pure lift for right handed players and φ = 90 o means pure slice).

Table 1 :

 1 Probability for the server to win the point as a function of serve rank (first, second) and court surface.

	surface	serve	win lose
	CLAY	first serve	66.28 33.72
		second serve 52.24 47.76
	GRASS	first serve second serve 54.85 45.15 71.19 28.81	gender serve	win lose
	HARD	first serve	68.34 31.66	women first serve	62.85 37.15
		second serve 52.68 47.32		second serve 49.43 50.57
	INDOORS first serve	72.01 27.99	men	first serve	71.00 29.00
		second serve 53.03 46.97		second serve 54.18 45.82

Table 2 :

 2 Probability for the server to win the point as a function of serve rank (first, second) and gender.

this is the average of each probability to win on first serve over the surfaces in Table1.

the x-axis is oriented along the depth of the court, the y-axis is parallel to the net.


Rupa Jagannathan

Sanja Petrovic

Angela Mckenna
email: angela.mckenna@nuh.nhs.uk

Louise Newton
email: louise.newton@nuh.nhs.uk

A Fuzzy Non-linear Similarity Measure for Case-Based Reasoning Systems for Radiotherapy Treatment Planning

Keywords: Case-based Reasoning, Fuzzy Logic, Radiotherapy Treatment Planning, Decision-Support Systems

This paper presents a decision support system for treatment planning in brain cancer radiotherapy. The aim of a radiotherapy treatment plan is to apply radiation in a way that destroys tumour cells but minimizes the damage to healthy tissue and organs at risk. Treatment planning for brain cancer patients is a complex decision-making process that relies heavily on the subjective experience and expert domain knowledge of clinicians. We propose to capture this experience by using case-based reasoning. Central to the working of our case-based reasoning system is a novel similarity measure that takes into account the non-linear effect of the individual case attributes on the similarity measure. The similarity measure employs fuzzy sets. Experiments, which were carried out to evaluate the similarity measure using real brain cancer patient cases show promising results.

Introduction

Radiotherapy is a treatment method that uses ionizing radiation in cancer control. The aim of radiotherapy treatment planning (RTP) is to deliver a tumouricidal radiation dose over the tumour region while minimizing the radiation received by healthy tissue and critical organs in the vicinity of the tumour. To realize this goal, a detailed treatment plan is created for each patient that describes exactly how a patient will be irradiated to achieve the required dose distribution. Parameters of RTP include the number of beams, the gantry angle of beams (the angle by which the apparatus applying the radiation beam is moved around the patient's body) and wedges that shape the beam. Oncologists, together with medical physicists, use their subjective experience and expert clinical knowledge to generate treatment plans. The RTP decision support system, under development in collaboration with the Nottingham University Hospitals, City Hospital campus, uses case-based reasoning to capture this experience and aid oncologists in the computation of plan parameters.

Case-based reasoning (CBR) is based on the concept that the solution of a problem can be derived from the solutions of similar problems [START_REF] Kolodner | Case-Based Reasoning[END_REF]. The main advantage of CBR in RTP over numerical optimisation [START_REF] Oldham | A Case Study Comparing the Relative Benefit of Optimizing Beam Weights, Wedge Angles, Beam Orientations and Tomotherapy in Stereotactic Radiotherapy of the Brain[END_REF] and rule-based methods [START_REF] Haas | Hybrid Optimisation Technique for Radiotherapy Treatment Planning[END_REF] is its capability to utilise non-quantifiable or subjective knowledge and benefit from the successes, errors and failures of previous treatment plans. In addition, since treatment plan generation is based on existing plans, new plans can be generated quickly and in line with the institution's preferences. Although, CBR has been widely applied in healthcare systems in the diagnosis and treatment of disease [START_REF] Holt | Medical Applications in Case-Based Reasoning[END_REF][START_REF] Schmidt | Cased-Based Reasoning for Medical Knowledge-Based Systems[END_REF], the research done on CBR in radiotherapy treatment planning has been limited. Case based reasoning has been used in dose planning for prostate cancer by Song et al [START_REF] Song | A Case-Based Reasoning Approach to Dose Planning in Radiotherapy[END_REF] and Mishra et al [START_REF] Mishra | A Knowledge-light Nonlinear Case-based Reasoning Approach to Radiotherapy Planning[END_REF]. The CBR system Roentgen, designed by Berger [START_REF] Berger | Roentgen: Radiation Therapy and Case-Based Reasoning[END_REF], aids radiotherapy planning for thorax cancer and retrieves cases from its case archive based on the similarity between patient geometric descriptors. However, no implementation details or experimental evaluation of the method are furnished. The work presented in this paper builds up on Berger's idea of using the geometry of the patient as input to the CBR system and implements this concept using a fuzzy non-linear similarity measure.

The choice of the similarity measure, which calculates how similar cases in the case base are to the target case, is fundamental. The popular nearest neighbour method [START_REF] Cover | Nearest Neighbor Pattern Classification[END_REF] matches each attribute in the target case to its corresponding attribute in the archive case. The aggregate similarity is given by the weighted sum of the individual similarities. However, if the similarity values with respect to individual attributes are very different in terms of their distribution or range they can not be directly compared or summed up to give an accurate representation of the aggregate similarity between two cases. Another common assumption of many CBR systems is that the similarity with respect to a single case attribute influences the aggregate similarity linearly over its range. In reality, however, the attribute similarity between two cases often behaves in a non-linear fashion [START_REF] Wang | Nonlinear Nearest-Neighbour Matching and its Application in Legal Precedent Retrieval[END_REF][START_REF] Cheng | A Fuzzy Inference System for Similarity Assessment in Case-Based Reasoning Systems: An Application to Product Design[END_REF]. That is, a high similarity value with respect to a non-linear attribute indicates a high suitability of the case solution to the target case but also a very low similarity of non-linear attribute has a detrimental effect on the suitability of the case. To overcome these limitations, our CBR system uses a novel fuzzy non-linear similarity measure. The performance of the proposed similarity measure is evaluated using real brain cancer patient cases obtained from the Nottingham City Hospital.

The paper is organized as follows. Section 2 describes the architecture of the CBR system and the non-linear fuzzy similarity measure. Section 3 presents some results obtained by comparing the performance of the standard weighted nearest neighbour and the fuzzy non-linear similarity measure. Section 4 discusses future research directions and concludes our work.

Architecture of the CBR System

Each case in the developed CBR system contains a patient description and a suitable treatment plan detailing the RTP beam configuration. The patient description consists of spatial information regarding the location of the planning target volume containing the tumour (PTV) and the organs at risk (OAR). The OAR commonly include the spinal cord, the eyes, the lens, the optic nerve and chiasm, the brainstem and the pituitary gland. The cases are sorted according to the OAR present in the vicinity of the PTV. Since procuring clinical patient data is a slow and difficult process, the case base currently contains only 24 brain cancer patient cases. However, we are continuously adding more cases to our case base as they become available.

Case Attributes

Case attributes identified to be relevant for similarity between two patients describe the location of the tumour and the spatial relationship between the tumour and the OAR. These attributes determine the geometry of a patient and are computed using 3D coordinates representing the PTV and OAR structures, which are extracted from the CT DICOM [START_REF]DICOM -Digital Imaging and Communications in Medicine[END_REF] image header files created for each patient. The following geometric descriptors are used.

• Volume of PTV, denoted by V.

• Distance between tumour and OAR, denoted by E, is defined as the minimum edge-to-edge distance between the tumour and the OAR. • Angle between tumour and OAR, denoted by A, is defined as the angle between the line connecting the origin of the image patient coordinate system and the tumour centroid and the line connecting the origin and the OAR centroid.

A Fuzzy Similarity Measure

A good definition of a similarity measure between two cases is of crucial importance for the retrieval process, which should retrieve from the case base a case with a treatment plan suitable for the target patient. A suitable treatment plan is one whose plan parameters achieve the desired radiation dose distribution in the target case patient. After consultation with medical physicists, we have arrived at the assumption that patients with similar geometric descriptors have similar treatment plans. The CBR system Roentgen is based on this assumption as well [START_REF] Berger | Roentgen: Radiation Therapy and Case-Based Reasoning[END_REF]. In the retrieval process, the similarity with respect to the geometric descriptors is calculated between the target case and those cases in the case base that contain the same OAR. The values of the three geometric descriptors are normalized to take values from the interval [0, 1]. Let C T be the target case and C C be a case from the case-base. The distance l d between attribute value T v in target case C T and attribute value C v in case C C with respect to attribute l , where l = A, E and V, is calculated as follows:

C T l v v d - = . (1)
Then, the similarity between C T and C C with respect to attribute l is:

l l d s - = 1 . ( 2
)
The similarity values with respect to each attribute have to be combined into an aggregate similarity value between two cases. In the nearest neighbour method, the aggregate similarity is generally defined as the weighted sum of the individual attribute similarities. This method however does not take into account the variations in the distribution of the similarity values of different attributes. In order to obtain an idea of the distribution of attribute similarities, we calculate the similarity between each case and every other case (with the same OAR) in our case base in a leave-one out fashion considering one attribute a time. Figure 1 shows the similarity calculated between all possible pairs of patient cases in the case base (arranged in decreasing order of similarity) for each attribute. The data suggests a trend in the distribution of the attribute similarity values between cases. We can see that the distribution and range of the similarity with respect to attributes A and E is comparable. However, it is different for attribute V. Hence, numerical similarity values with respect to attribute A or E do not necessarily have the same meaning as an equal numerical similarity value with respect to attribute V. To illustrate this point, consider the case pair C 10 and C 23 .

In figure 1, we can see that the similarity between C 10 and C 23 with respect to the volume V is '0.72', which in linguistic terms can be considered as an 'average' volume similarity. However, if the similarity between C 10 and C 23 with respect to attribute angle was '0.72', we observe that this value is considered as a relatively low angle similarity compared to all other similarity values for that attribute and we conclude that the two cases are not very similar with respect to the attribute angle. That is, the same numerical similarity value has different meanings for different attributes. Thus, we can not compare directly the numerical values of the similarity and a simple sum of individual attribute similarities does not necessarily give an accurate representation of the similarity between two cases. In order to sensibly compare attribute similarities and generate an aggregate similarity measure, the CBR system has to understand what numerical similarity value actually constitutes a "high", "average" or "low" similarity for each attribute. We therefore propose the use of fuzzy sets. We define the fuzzy sets High, Avg and Low, which denote low similarity, average similarity and high similarity, respectively, for each attribute. The fuzzy membership functions of the three sets are defined for each attribute based on the minimum, maximum and average of the corresponding similarity values found across the case base. They, therefore, give a realistic indication of what constitutes a relatively "high similarity", "average similarity" or "low similarity" for an attribute.

= = Μ 3 V , E , A l s , l l s w µ . ( 3
)
where l w denotes the weight of attribute l , l = A, E, V and l,s , s = High, Avg, Low, is the membership degree of the attribute similarity to the fuzzy sets High, Avg and Low.

A large value of component

M w M w M w s - - = . (4)
The weights w High , w Avg and w Low determine the importance of the fuzzy sets High, Avg and Low respectively.

Similarity

M embership

Low Avg High

Fig. 2. Fuzzy membership functions for sets High, Low and Avg.

Non-linear Similarity

Preliminary experiments have shown that the influence of the Angle similarity (s A ) on the aggregate similarity is not linear. That is, if s A is high, then it is very likely that the treatment plan of the retrieved case is very suitable and contains a beam configuration that would achieve the desired radiation dose distribution for the target patient. However, if s A is low, then the corresponding treatment plan is very likely to be rather unsuitable for the target case. In other words, to retrieve cases whose treatment plans are very suitable for the target case, the aggregate similarity between cases has to be increased if the value s A is high and reduced or penalised if s A is low. This effect is not covered by the standard method of multiplying weights with the similarity value assigned to attributes, since a high weight for attribute A, would increase the aggregate similarity value not only if s A is high but also if s A is low.

For the similarity with respect to V, s V , the effect is non-linear too, but in an opposite fashion. If s V is high, it does not substantially increase the likelihood that the treatment plan of the retrieved case is suitable. Similarly, if s V is low, it does not necessarily mean that the corresponding treatment plan is highly unsuitable i.e. the influence of both high and low similarity due to attribute V is reduced or smoothed.

Expression 5 presents the aggregate similarity measure, which takes into account the non-linear effect of attributes, where l k is the non-linearity factor of attributel. We can see that the fuzzy similarity measure, wFuzzy, has a lower average error (18.6°) than the nearest neighbour method, wNN (15.5°). The error is further reduced with the fuzzy non-linear similarity measure, wNL (14.7°). Following consultation with hospital staff, an angle error of 20° per beam is acceptable for most cases. Figure 3b, shows the number of cases for which the similarity measure retrieves a case, in which the average error per beam is less than 20°. We can see that the percentage of cases within this beam angle tolerance is 58% for wNN, 75% for wFuzzy and 79% for wNL similarity measure.

Conclusion

In this paper, we described a case-based reasoning system (CBR) that is under development to aid oncologists with radiotherapy treatment planning in brain cancer.

The results confirm the assumption that similar treatment plans can be derived from similar geometric patient descriptors, as long as these descriptors are combined in an appropriate fashion in the similarity measure. Therefore, a key characteristic of our system is the fuzzy non-linear similarity measure that computes the similarity between cases in order to retrieve the case with the most suitable treatment plan. The fuzzy nature of the similarity measure allows aggregation of the similarity values with respect to single attributes to provide an accurate representation of the total similarity while taking into account the non-linear behaviour of attributes. The method was compared to the commonly used nearest neighbour method. The experimental results using real brain cancer patient cases show an improvement with the fuzzy non-linear similarity measure. The attribute weights and the non-linearity factor k have been set empirically, but the performance of the similarity measure can be improved by optimising the weights. A detailed weight analysis to identify the optimum values of the attribute weights and k is currently in process. We are aiming to find a set of weights that return the lowest average error with the leave-one-out cross validation method. Preliminary results show a significant reduction in beam angle error with optimized attribute weights.

In CBR systems the performance of the case retrieval mechanism and its similarity measure is crucial to the working of a CBR system. However, frequently, differences do exist between the target case and the retrieved case. Once the most similar case has been retrieved, its solution generally has to be adapted to the specific needs of the target case. Adaptation can be done by adjusting the beam configuration according to the geometric displacement in the location of the tumour and OAR structures of the target case compared to the retrieved case. Another method evaluates the plan using a treatment planning system to identify dose violations of OAR and confirm tumour coverage. The adaptation module of the CBR system then has to tweak the plan parameters of the retrieved case to resolve the violations. Adaptation can be carried out using if-then rules or by using another case-based reasoning inference mechanism.

Fig. 1 .

 1 Fig. 1. The graph shows the similarity values obtained between all possible pairs of cases.

Figure 2

 2 Figure 2 presents the fuzzy membership functions defined for attributes A, E and V. The membership degree s , l µ of similarity s l with respect to attribute l is computed for each fuzzy set s, High, Avg and Low. The aggregate similarity consists of the High, Avg and Low component s Μ defined as the sum of the membership degrees of the attribute similarities to the fuzzy sets High, Avg and Low of the corresponding attributes, shown in expression 3.

Fig. 3 .

 3 Fig. 3. a) The average error between the beam angles of the treatment plan of the retrieved case and the actual treatment plan. b) The percentage of cases, for which the similarity measure retrieves a treatment plan with acceptable beam angle difference.

Acknowledgements. The authors would like to thank the Engineering and Physics Science Research Council (EPSRC), UK, (Ref No. EP/C549511/1) and the City Hospital, Nottingham University Hospitals NHS Trust (UK) for their support.

(

) ( )

To achieve a stretching non-linear effect, k should be smaller than 1, to achieve a smoothing non-linear effect, k should be greater than 1 and k should be equal to 1 if no non-linear effect is desired. After initial experiments, we have used

, for attribute E.

Experimental Results

We tested the performance of the fuzzy non-linear similarity measure using brain cancer patient cases. Since all cases contain the treatment plan used in that case, we can evaluate the similarity measure by comparing the treatment plan of the retrieved case with the actual treatment plan. In the experiments, we looked mainly at the beam configuration of the plan, in particular the beam angles. The beam angles consider the position of the gantry, the patient couch and the beams. The difference in the beam angle values in the plan of the retrieved case and the actual treatment plan of the target case gives the error. The error is averaged over all cases to give an indication of the quality of the similarity measure. A commonly employed evaluation strategy in CBR is the leave-one out cross validation method. Each case in the case base is consecutively made the target case and the most similar case among the remaining cases in the case base is retrieved. In the future when more cases are available, the case base can also be divided into training and test cases. Figure 3a shows the average errors in angle degrees per beam across all cases obtained when using the weighted nearest neighbour similarity measure (wNN), the fuzzy similarity measure (wFuzzy) given in expression (4) and the fuzzy non-linear similarity measure (wNL) given in expression [START_REF] Schmidt | Cased-Based Reasoning for Medical Knowledge-Based Systems[END_REF]. Following discussion with hospital staff, attribute A is the most important, followed by E and V. Therefore, the attribute weights are set at w A = 0.6, w E = 0.3 and w V = 0.1 for all similarity measures. The weights for the fuzzy sets High, Low and Avg are set at wHigh =3, wLow = 2 and wAvg =1.