Shanyu Tang
email: s.tang@londonmet.ac.uk

Hassan B Kazemian

Simulation of web data traffic patterns using fractal statistical modelling

Keywords: data traffic, patterns, fractal modelling

This paper describes statistical analysis of web data traffic to identify the probability distribution best fitting the connection arrivals and to evaluate whether data traffic traces follow heavy-tail distributions -an indication of fractal behaviour, which is in contrast to conventional data traffic.

Modelling of the fractal nature of web data traffic is used to specify classes of fractal computing methods which are capable of accurately describing the burstiness behaviour of the measured data, thereby establishing web data traffic patterns.

Introduction

Data traffic is a vibrant area of queuing research. When modelling network traffic, network arrivals such as packet and connection arrivals are often modelled as Poisson processes for analytical simplicity because such processes have attractive theoretical properties [START_REF] Frost | Traffic modeling for telecommunications networks[END_REF]. However, a number of studies have shown that for both local-area [START_REF] Gusella | A measurement study of diskless workstation traffic on an Ethernet[END_REF][START_REF] Fowler | Local area network traffic characteristics, with implications for broadband network congestion management[END_REF] and wide-area [START_REF] Danzig | An artificial workload model of TCP/IP internetworks[END_REF][START_REF] Paxson | Wide area traffic: the failure of Poisson modeling[END_REF] network traffic, packet interarrivals are clearly not exponentially distributed as expected from the Poisson processes. The distinctly non-Poisson nature of the arrival process has led to a substantial effort in developing statistical models of the traffic. These statistical models differ from the renewal and Markov-chain-type models that formed the basis for traditional queuing models [START_REF] Heyamn | Some issues in performance modeling of data traffic[END_REF], so they inspired new queuing models and solution methods. Previous wide-area traffic studies were largely focused on FTP, TELNET, NNTP and SMTP (email) traffic [START_REF] Danzig | An artificial workload model of TCP/IP internetworks[END_REF][START_REF] Paxson | Wide area traffic: the failure of Poisson modeling[END_REF][START_REF] Heyamn | Some issues in performance modeling of data traffic[END_REF][START_REF] Paxson | Empirically-Derived Analytic Models of Wide-Area TCP Connections[END_REF], with little attention being given to web traffic [START_REF] Nabe | Analysis and modeling of world wide web traffic for capacity dimensioning of internet access lines[END_REF]. The aim of this work is to tackle web traffic by providing a view of web data patterns. Since web traffic accounts for a large proportion of the traffic on the Internet, understanding the nature of web traffic is increasingly important. Like other wide-area traffic, web traffic may also have burstiness behaviour; i.e. web traffic is burst on many or all time scales. This burstiness behaviour is analogous to the self-similar or fractal-like behaviour, which exists in many natural phenomena such as Brownian motion, turbulent flow, atmospheric pressure, the distribution of stars and the activity of the stock market [START_REF] Feder | Fractals[END_REF][START_REF] Falconer | Fractal Geometry: Mathematical Foundations and Applications[END_REF][START_REF] Kaye | A Random Walk Through Fractal Dimensions[END_REF][START_REF] Tang | Computer simulation of fractal structure of flocs[END_REF][START_REF] Rezaul | Empirical distribution function test for World Wide Web traffic[END_REF], which are much better characterised by fractal geometry theory than by Euclidean geometry.

Performance evaluation is important for assessing the effectiveness of data traffic pattern prediction, besides monitoring and verifying compliance with network performance goals [START_REF] Awduche | RFC 3272: Overview and Principles of Internet Traffic Engineering[END_REF]. Results from performance evaluation are used to identify existing problems, guide network re-optimisation, and aid in the prediction of potential future problems. Performance evaluation can be conducted in many different ways. The most notable techniques include analytical methods, simulation, and empirical methods based on measurements.

Previous research found that the light-and heavy-tailedness of web data traffic patterns. The tail weight was estimated by using the Hill method in most of the cases presented in [START_REF] Judith | From Network Management Collection to Traffic performance modelling: Challenges and lessons learned[END_REF][START_REF] Crovella | Heavy-Tailed Probability Distributions in the World Wide Web[END_REF][START_REF] Labs | Busier Networks create smoother traffic flow[END_REF]. The method is a statistical technique for analysing heavy-tailed phenomena. The less heavy-tailed the data traffic, the more benefit of a buffer [START_REF] Clark | Provisioning for Bursty Internet Traffic: Implications for industry and Internet structure[END_REF]. Adding multimedia files to the set of text files led to the increase in the weight of the tail and the distribution of text files might itself be heavy-tailed [START_REF] Crovella | Heavy-Tailed Probability Distributions in the World Wide Web[END_REF]. For light-tailed input, the delay distribution has an exponential tail; whereas the delay distribution follows a lognormal distribution for heavy-tailed input [START_REF] Squillante | Internet Traffic: Periodicity, Tail Behaviour, and Performance Implications[END_REF]. A hypothetical example [START_REF] Prophet Statguide | Goodness of Fit (Chi-square) Test[END_REF] of a normal probability plot for data sampled from a distribution has been illustrated from where the light-tailedness and the heavytailedness of the data traffic pattern can be identified by visualization.

Characterising heavy-tail behaviour of web data traffic

Several probability distributions including conventional exponential and normal distributions have been chosen to fit the web traffic traces e.g. web document size, request interval and access frequencies collected at the measurement points. Previous work [START_REF] Paxson | Wide area traffic: the failure of Poisson modeling[END_REF][START_REF] Heyamn | Some issues in performance modeling of data traffic[END_REF] shows that web traffic such as telnet and ftp followed heavytail distributions. This means that values are distributed over a very wide range and that the larger values, even if less probable, may still account for a significant portion of the traffic. Several long-tailed distributions that were chosen to fit the traces are a lognormal distribution, an extreme distribution, a log-extreme distribution, a Weibull distribution and a Pareto distribution.

Empirical distribution function (EDF) test was also used to decide whether the data (traces) were from a particular distribution and best fittings were made using the EDF statistics. Parameters (e.g. location parameter, scale parameter, shape parameter) that affect in fitting the curve are highlighted to realize their effects. The analytical expressions of some distributions are detailed in the following sections.

Log-normal distribution

The lognormal law or log-Gaussian frequency function is defined as x is the size corresponding to the 16% cumulative frequency, and 84

          g g g x x x d d     2 2 log
x is the size corresponding to the 84% cumulative frequency.

Pareto distribution

The cumulative distribution for the Pareto (simple) random variable is obtained from

                 x x x X P x F , 0 , ; 1 ) ( ) (
where  is the location parameter, and  is the shape parameter.

EDF test

The EDF is a step function calculated from the sample, estimating the population distribution function. EDF statistics are measures of the discrepancy between the EDF and a given distribution function, and they are used for testing the fitting of the sample to the distribution. Suppose a given random sample of size n is X 1 , …, X n and X (1) < X (2) < … < X (n) be the order statistics; suppose further that the distribution of x is F(x). The empirical distribution function, F n (x), is given by

         x n x ns observatio of number x F n ;
The calculation can be conducted by using the Probability Integral Transformation (PIT), Z = F(X); when F(x) is the true distribution of x, the new random variable Z is uniformly distributed between 0 and 1. Then Z has the distribution function F*(z) = z, 0  z  1. The following formulas can be used for calculating EDF statistics from the Z-values. The formulas involve the Z-values arranged in an ascending order, Z (1) <Z (2) < … < Z (n) .

                                                        i n i n i n i i n i i Z Z i n n A n Z mean sample z where z n W U n n i Z W 1 1 2 1 2 2 2 2 1 2 1 log log 1 2 1 ; 5 . 0 ) 12 /( 1 )} 2 /( ) 1 2 ( {
where log x represents log e x (natural logarithm).

The decision can be made by comparing the calculated value with the tabulated value [START_REF] Fowler | Local area network traffic characteristics, with implications for broadband network congestion management[END_REF]. The hypothesis is rejected at significance level p if the calculated value is greater than the tabulated value given for level p [START_REF] D'agostino | Tests based on EDF statistics[END_REF].

Data traffic collection

The first step in understanding network traffic is the collection of trace data. The LBL-CONN-7 trace [22] was collected at the Lawrence Berkeley Laboratory (LBL), located in Berkeley, California, containing thirty days' worth of all wide-area TCP connections between the LBL and the rest of the world. The reduced trace was generated by TCP-reduce. TCP-reduce is a collection of Bourne shell scripts for reducing tcpdump traces to one-line summaries of each TCP connection present in the trace. The scripts are TCP-reduce, which takes a tcpdump trace file as an argument and writes a sorted summary to stdout, TCP-conn (an internal awk script that does all the work) and TCP-summary (an awk script that generates a per-protocol summary of all the TCP connections produced by TCP-reduce). The scripts were written using Bourne shell, tcpdump and the common Unix utilities sed, sort and awk. The trace was written as an ASCII file with one line per connection with the columns such as timestamp, duration, protocol, bytes sent by originator of the connection, bytes sent by responder to the connection, local host, remote host, state that the connection ended in and flags.

The trace ran from midnight, Thursday, 16 September 1993 through midnight, Friday, 15 October 1993 (times are Pacific Standard Time), capturing 606,497 wide-area connections. The tracing was performed on the Ethernet DMZ network over which flows all traffic into or out of the LBL. The raw trace was made using tcpdump on a Sun SPARC station using the BPF kernel packet filter. Fewer than 15 SYN/FIN/RST packets in a million were dropped. Timestamps had microsecond precision. The traffic was filtered to exclude connections with nearby UCB except for nntp.

A special care was taken to check the sanity of the data, as any irregularity or mistakes could set back the entire analysis for an extended period of time. The original data file was divided into smaller file and then the size (byte) of the file was arrayed into different bins for convenience of the compactness of the data file. It is very congenial to array the file size into a bin as thousands of samples can be analysed in an expected range.

Results and Discussion

Traffic data in byte sent by originator

Fig. 1(a) shows the web file size distribution in percentage. The highest percentage of the file size transferred is 23.4% at 500 bytes whereas the lowest is 0.30% at 23,000 bytes. The lognormal distribution does not follow the trace data satisfactorily. The t-test for one sample (web file size) in two tails distribution found to be 3.4916, which is greater than the tabulated value (Table 1), indicating that the trace data distribution may be tailed. Fig. 1(b) illustrates comparisons of the cumulative distribution between the web file data and the data generated using Pareto model. The model data somewhat fit the web data within the range 10 to 1000 byte, which is an indication of burstiness behaviour. The test statistics (A2, W2 and U2 tests) confirm (Table 2) that the trace data follow the normal distribution. There is contradiction between EDF test statistics and t-test as the t-test does not support the trace data for normality. The A2 test rejected the hypothesis for Weibull distribution whereas W2 and U2 tests show that the data are significant at the level  = 0.025 which is consistent with the results shown in Fig. 2(c). The A2 test also rejected the hypothesis for generalized Pareto distribution where the data are significant at the level of 0.001, whereas W2 accepted the hypothesis at the level of 0.025 that makes sense with the graphical results shown in Fig. 2(b). The three test statistics show that the trace data are significant at  = 0.10 and =0.25 for extreme value and exponential distributions respectively, also  = 0.05 (shown by A2 and W2 tests) and  = 0.10 (by U2 test) for logistic distribution. The three distributions could have shown better fitness than others according to significance level. Unfortunately, poor matches were observed as shown in Figs. 2(a),(d) and (e) which is questionable. Fig. 2(f

Traffic data in byte sent by responder

Fig. 3(a) shows the web file size distribution in percentage. The highest percentage of the file size transferred is 19.48% at 150 bytes whereas the lowest is 0.12% at 600 bytes. The lognormal distribution only partly fits the web data. The t-test for one sample (web file size) in two tails distribution found to be 7.3395, which is greater than the tabulated value (Tables 1), indicating that the trace data do not follow normal distribution. The test statistics (A2, W2 and U2 tests) confirm (Table 3) that the trace data follow the normal distribution. The test statistics also show that the data are significant at the levels of 0.10, 0.15, 0.25, and 0.05 for extreme value, Weibull, logistic and Pareto (generalized) distributions, respectively, which is consistent with the results shown in Figs. 4(b), (c), (d), and (e). For exponential distribution, the A2 and W2 tests show the data are significant at the level  = 0.05, and U2 shows at the level of 0.10 and the effect of these significance levels observed in Fig. 4(a) which is understandable. Likewise, there is contradiction between EDF test statistics and t-test as t-test does not support the trace data for normality.

Web File Size Distribution

Conclusions

The T-test and EDF test statistics (A2, W2 and U2) have been used to evaluate how well a particular distribution describes the real web traffic data. The models employed in the work include lognormal, Pareto (simple), exponential, Pareto (generalized), Weibull (three parameters), logistic and extreme value distributions.

There is contradiction between EDF test statistics and t-test; the t-test results do not support the web trace data for normality, but the EDF test accepts the hypothesis that the trace data follow normal distribution. The discrepancy may be due to the different sampling methods, which has to be investigated in the future.

There are satisfactory fittings of curves observed between Weibull and the generalized Pareto model data and the real trace data, which has been confirmed by the EDF test for the two distributions. The results show that Weibull (three parameters) is the most suitable model to approximate the web traffic. In addition, the generalized Pareto model is more suitable for analysing traffic fractal-like behaviour (burstiness behaviour) than simple Pareto model.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Web file size distribution and Cumulative distribution of LBL_CONN_7 (byte sent by originator) data

  ) illustrates the combination of Figs. 2(a), (b), (c), (d), and (e) from where the best fitness can be deemed by comparison.

Fig. 3 .Fig. 4 Fig. 4 .

 344 Fig. 3. Web file size distribution and Cumulative distribution of LBL_CONN_7 (byte sent by responder) data

Fig. 4 (

 4 Fig. 4(f) shows the combination of Figs 4(a), (b), (c), (d), and (e) from where comparisons of goodness-of-fit of trace data to different distributions can be made.

Table 1 .

 1 T-test results for web file size data

	Traffic	t	Degrees of	t0.05
		(Experimental)	freedom	(tabulated value)
	LBL-CONN-7-	3.4916	11	2.201
	originator
	LBL-CONN-7-	7.3395	12	2.179
	responder

Acknowledgements

The authors would like to thank LBL for providing the web traffic data.


Manuel Cruz

Jorge M Santos

Nuno Cruz
email: nbdfcruz@gmail.com

Maximum Shear Modulus Prediction by Marchetti Dilatometer Test Using Neural Networks

The use of Neural Networks for modeling systems has been widespread, in particular within areas where the great amount of available data and the complexity of the systems keeps the problem very unfriendly to treat following traditional data analysis methodologies. In the last two decades, small strain shear modulus became one of the most important geotechnical parameters to characterize soil stiffness. Finite element analysis have shown that in-situ stiffness of soils and rocks is much higher than was previously thought, and that stress-strain behaviour of these materials is non-linear in most cases with small strain levels, especially in the ground around retaining walls, foundations and tunnels typically in the order of 10 -2 to 10 -4 of strain. Although the best approach seems to be based in measuring seismic wave velocities, deriving the parameter through correlations with in-situ tests is usually considered very useful for design practice. In this work, a new approach using Neural Networks is proposed for sedimentary soils and the results are discussed and compared with some of the most common available methodologies for this evaluation.

Introduction

Maximum shear modulus, G 0 , has been increasingly introduced in stiffness evaluations for design purposes, for the last 10-15 years. At the present moment, one of the best ways for measuring it is to evaluate compression and shear wave velocities and thus obtain results supported by theoretical interpretations. The advantages of this approach are widely known, mainly because test readings are taken through intact soil in its in-situ stress and saturation levels, thus practically undisturbed, and also dynamic stiffness can be close to operational static values [START_REF] Clayton | Stiffness of geomaterials at very small strains[END_REF][START_REF] Fahey | Soil stiffness values for foundation settlement analysis[END_REF]. However, the use of seismic measures implies a specific test and so, many atempts have been made to correlate other in-situ test parameters, such as those obtained by Standard Penetration Test, SPT [START_REF] Peck | Foundation Engineering[END_REF], Piezocone Test, CPTu [START_REF] Lunne | Cone penetration testing in geotechnical practice[END_REF] or Marchetti Dilatometer Test, DMT [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF][START_REF] Cruz | Dmt experience in iberian transported soils[END_REF][START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF] with G0, using cross calibrations (e.g.,seismic or triaxial among others).

The DMT test is one of the most appropriated for this task (although some relations can be settled using SPT or CPTu) since it uses a measurement of a load range related with a specific displacement, which can be used to deduce highly accurate stress-strain relationship (E D ), supported by Theory of Elasticity. Moreover, the type of soil can be numerically represented by DMT Material Index, I D , while in situ density, overconsolidation ratio (OCR) and cementation influences can be represented by lateral stress index, K D , allowing for high quality calibration of the basic stress-strain relationship [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF].

The present paper is organized as follows: In Section 2 an overview of the most important correlations used in this subject is reported, in order to get a clear view of the problem context. In Section 3 the data set is presented, as well as the methodology used for the Neural Networks application on this subject and its results and subsequent discussion. As usual, the last section will be used for conclusions and final remarks.

G 0 prediction by DMT

Marchetti dilatometer test or flat dilatometer, commonly designated by DMT, has been increasingly used and it is one of the most versatile tools for soil characterization, namely loose to medium compacted granular soils and soft to medium clays, or even stiffer if a good reaction system is provided. The test was developed by Silvano Marchetti [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF] and can be seen as a combination of both Piezocone and Pressuremeter tests with some details that really makes it a very interesting test available for modern geotechnical characterization [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF]. The main reasons for its usefulness deriving geotechnical parameters are related to the simplicity (no need of skilled operators) and the speed of execution (testing a 10 m deep profile takes around 1 hour to complete) generating quasi-continuous data profiles of high accuracy and reproducibility.

In its essence, dilatometer is a stainless steel flat blade (14 mm thick, 95 mm wide and 220 mm length) with a flexible steel membrane (60 mm in diameter) in one of its faces. The blade is connected to a control unit on the ground surface by a pneumatic-electrical cable that goes inside the position rods, ensuring electric continuity and the transmission of the gas pressure required to expand the membrane. The gas is supplied by a connected tank/bottle and flows through the pneumatic cable to the control unit equipped with a pressure regulator, pressure gauges, an audio-visual signal and vent valves. The equipment is pushed (most preferable) or driven into the ground, by means of a CPTu rig or similar, and the expansion test is performed every 20cm. The (basic) pressures required for liftoff the diaphragm (P 0 ), to deflect 1.1mm the centre of the membrane (P 1 ) and at which the diaphragm returns to its initial position (P 2 or closing pressure) are recorded. Due to the balance of zero pressure measurement method (null method), DMT readings are highly accurate even in extremely soft soils, and at the same time the blade is robust enough to penetrate soft rock or gravel (in the latter, pressure readings are not possible), supporting safely 250kN of pushing force. The test is found especially suitable for sands, silts and clays where the grains are smaller (typically 1 10 to 1 5 ) compared to the membrane dimension [START_REF] Marchetti | The flat dilatometer: Design applications[END_REF]. Four intermediate parameters, Material Index (I D ), Dilatometer Modulus (E D ), Horizontal Stress Index (K D ) and Pore Pressure Index (U D ), are deduced from the basic pressures P 0 , P 1 and P 2 , having some recognizable physical meaning and some engineering usefulness [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF], as it will be discussed below. The deduction of current geotechnical soil parameters is obtained from these intermediate parameters covering a wide range of possibilities. In the context of the present work only E D , I D and K D have a physical meaning on the determination of G 0 , so they will be succinctly described as follows [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF]:

1. Material Index, I D : Marchetti [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF] defined Material Index, I D , as the difference between P 1 and P 0 basic pressures measured normalized in terms of the effective lift-off pressure. The I D parameter is one of the most valuable indexes deduced from DMT, due to its ability to identify soils throughout a numerical value that can be easily introduced in specific formulae for deriving geotechnical parameters. In a simple form, it could be said that I D is a "fine-content-influence meter" [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF], providing the interesting possibility of defining dominant behaviours in mixed soils, usually very difficult to interpret when only grain size is available, thus it may be associate to an index reflecting an engineering behaviour. 2. Horizontal Stress Index, K D : The horizontal stress index [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF] was defined to be comparable to the at rest earth pressure coefficient, K 0 , and thus its determination is obtained by the effective lift-off pressure (P 0 ) normalized by the in-situ effective vertical stress. K D is a very versatile parameter since it provides the basis to assess several soil parameters such as those related with state of stress, stress history and strength, and shows dependency on several factors namely cementation and ageing, relative density, stress cycles and natural overconsolidation resulting from superficial removal, among others. The parameter can be regarded as a K 0 amplified by penetration effects [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF] and displays a typical profile very similar in shape to OCR profiles, giving useful information not only about stress history but also on the presence of cementation structures [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF]. Since undrained shear strength of fine soils can be related and obtained via OCR and the relation between K 0 and angle of shearing resistance is well stated by soil mechanics theories, then the parameter is also used with success in deriving shear strength. 3. Dilatometer Modulus, E D : Stiffness behaviour of soils is generally represented by soil moduli, and thus the base for in-situ data reduction. Stiffness behaviour of soils is generally represented by soil moduli, and thus the base for in-situ data reduction. Theory of Elasticity is used to derive dilatometer modulus, E D [START_REF] Marchetti | In-situ tests by flat dilatometer[END_REF] , by considering that membrane expansion into the surrounding soil can be associated to the loading of a flexible circular area of an elastic half-space, and thus the outward movement of the membrane centre under a normal pressure (P 1 -P 0 ) can be calculated. In short, E D is a parameter that includes both Young modulus (E) and Poisson's coefficient (ν)

and can be expressed as follows:

E D = E 1 -ν 2 = 34.7(P 1 -P 0 ) (1)
Generally speaking, soil moduli depend on stress history, stress and strain levels drainage conditions and stress paths. The more commonly used moduli are constrained modulus (M ), drained and undrained compressive Young modulus (E 0 and E u ) and small-strain shear modulus (G 0 ), this latter being assumed as purely elastic and associated to dynamic low energy loading.

Maximum shear modulus, G 0 , is indicated by several investigators [START_REF] Fahey | Soil stiffness values for foundation settlement analysis[END_REF][START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF][START_REF] Mayne | Interrelationships of dmt and cpt in soft clays[END_REF] as the fundamental parameter of the ground. If properly normalized, with respect to void ratio and effective stress, could be seen as independent of the type of loading, number of loading cycles, strain rate and stress/strain history [START_REF] Mayne | Interrelationships of dmt and cpt in soft clays[END_REF]. It can be accurately deduced through shear wave velocities,

G 0 = ρv 2 s (2)
where ρ stands for density and v s for shear wave velocity.

However, the use of a specific seismic test imply an extra cost, since it can only supply this geotechnical parameter, leaving strength and insitu state of stress information dependent on other tests. Therefore, several attempts to model the maximum shear modulus as a function of DMT intermediate parameters for sedimentary soils have been made in the last decade. Hryciw [START_REF] Hryciw | Small-strain-shear modulus of soil by dilatometer[END_REF] proposed a methodology for all types of sedimentary soils, developed from indirect method of Hardin & Blandford [START_REF] Hardin | Elasticity of particulate materials[END_REF]. Despite its theoretical base, references on this method are scarce and in general it is not applied for practical analysis due to some scatter around the determination, as illustrated by the results obtained in Portuguese normally consolidated clays [START_REF] Cruz | Dmt experience in iberian transported soils[END_REF]. The reasons for this scatter may be related to an expectable deviation of K 0 due to the important disturbance effects generated by penetration. On the other hand, this methodology ignores dilatometer modulus, E D , commonly recognized as a highly accurate stress-strain evaluation, and also lateral stress index, K D , and material index, I D , which are the main reasons for the accuracy in stiffness evaluation offered by DMT tests [START_REF] Cruz | Dmt experience in iberian transported soils[END_REF]. Being so, the most common approaches [START_REF] Jamiolkowski | New developments in field and laboratory testing of soilsladd, c.c[END_REF][START_REF] Sully | Correlation of maximum shear modulus with dmt test results in sand[END_REF][START_REF] Tanaka | Characterization of sandy soils using CPT and DMT[END_REF] with reasonable results concentrated in correlating directly G 0 with E D or M DM T (constrained modulus), which have revealed linear correlations with slopes controlled by the type of soil. In 2006, Cruz [START_REF] Cruz | Dmt experience in iberian transported soils[END_REF] proposed a generalization of this approach, trying to model the ratio

R G ≡ G0
ED as a function of I D . In 2008, Marchetti [START_REF] Marchetti | From Research to Practice in Geotechnical EngineeringD[END_REF] using the commonly accepted fact that maximum shear modulus is influenced by initial density and considering that this is well represented by K D , studied the evolution of both R G and G 0 /M DM T with K D and found different but parallel trends as function of type of soil (that is I D ), recommending the second ratio to be used in deriving G 0 from DMT, as consequence of a lower scatter. In 2010, using the Theory of Elasticity, Cruz [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF] approximate G 0 as a non-linear function of I D , E D and K D , from where a promising median of relative errors close to 0.21 with a mean(std) around 0.29(0.28) were obtained. It is worth mention that comparing with the previous approach -R G -this approximation, using the same data, lowered the mean and median of relative errors in more than 0.05 maintaining the standard deviation (Table 2).

In this work, to infer about the results quality it will be used some of the same indicators used by Hryciw, Cruz and others that are: the median, the arithmetic mean and standard deviation of the relative errors

δ i G0 = | G 0 (i) -G 0 (i)| |G 0 (i)| ; i = 1, 2, ..., N (3)
where G 0 (i) stands for the predicted value and G 0 (i) for the measured value seismic wave velocities (which is assumed to be correct). A final remark to point out that since in this work the no-intercept regression is sometimes used, the R 2 values will not be presented as they can been meaningfull in this case [START_REF] Huang | Draper Transformations, regression geometry and R2[END_REF]. It is also worth to remark that in the context of DMT and from the engineering point of view, median is the parameter of choice for assessing the model quality [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF] since the final value for maximum shear modulus relies on all set of results obtained in each geotechnical unit or layer.

Data Set, Experiments and Results

The data set used in the model is exactly the same used in the development of the non-linear G 0 approximation used by Cruz [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF], resulting from 860 DMT measurements performed in Portugal by Cruz and world wide by Marchetti et al. [START_REF] Marchetti | From Research to Practice in Geotechnical EngineeringD[END_REF] (data kindly granted by Marchetti for this purpose), which included data obtained in all kinds of sedimentary soils, namely clays, silty clays, clayey silts, silts, sandy silts, silty sands and sands. The main statistical measures of I D , E D , K D and G 0 parameters are given in Table 1 as well as in Figure 1. Several types of neural networks (NN) were used in order to improve the results obtained with traditional approaches and to achieve the best results with this kind of tool. The purpose was to find the best NN algorithm for the DMT problem. This process started by performing some experiments with traditional Multi Layer Perceptrons (MLP's) with different learning algorithms, namely Quasinewton, Conjugated Gradient and Scaled Conjugated Gradient (SCG) but, as it will be shown latter, the results were not as good as expected, and so Radial Basis Function (RBF), Bayesian and Fitting neural networks were then applied. The results obtained with this neural networks were slightly better than those with MLP's but still not significantly better. The last approach was to use Support Vector Regression in order to improve the results with the previous NN. For all the experiments the 10 fold cross validation method with 20 repetitions was used, since this is the most common and widely accepted methodology to guarantee a good neural network generalization. For each NN a huge set of experiments was performed, varying the involved parameters such as the number of neurons in the MLP hidden layer, the number of epochs or the minimum error for stopping criteria. The results presented in Tables 3 to 7 are therefore the best ones for each regression algorithm and represents the mean of the 10 × 20 performed tests for each best configuration. It is important to stress the fact that, when compared to traditional approaches where all the data is used to build the model, this methodology tends to produce higher standard deviations since in each experiment only a fraction of the available data is used to evaluate the model.

Values
As stated in Section 2, Cruz [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF] was able to improve older results applying non-linear function approximation. In Table 2 the median and mean (standard deviation) of the relative errors obtained by that approaches are presented.

Table 3 shows the best results for the first performed experiments with common MLP regression using three different learning algorithms: Quasinewton, Conjugated Gradient and SCG. It can be seen that there was a small improvement in the results when compared to those of Table 2 but not very significant.

Type

Median/Mean(std)

Non-Linear Regression G0 = α ED (ID) β 0.28/0.34(0.29) G0 = ED + ED e (α + βI D + γ log(K D )) 0.21/0.29(0.28) The second set of experiments was performed using Bayesian, RBF's (with thin plate spline activation function) and Fitting neural networks regression, by applying several combinations for the number of iterations in RBF's or the number of inner and outer loops for the Bayesian NN's. The best obtained results with these algorithms are presented in Table 4. These results are not significantly better than those obtained with common MLP's, with the exception of the fitting neural network that shows a significant improvement when compared with the others and also with traditional approaches.

An attempt to improve the quality of results was carried out by using Support Vector Regression (SVR). Support Vector Machines [START_REF] Cortes | Support-vector Networks[END_REF] are based on the statistical learning theory from Vapnik and are specially suited for classification. However, there are also algorithms based in the same approach for regression problems known as Support Vector Regression (SVR). Two different kinds of SVR algorithms: -SVR, from Vapnik [START_REF] Vapnik | Statistical Learning Theory[END_REF] and ν-SVR from Schölkopf [START_REF] Schölkopf | New support vector algorithms[END_REF] were applied, which differ in the fact that ν-SVR uses an extra parameter ν ∈ (0, 1] to control the number of support vectors. For these experiments different values for the parameter C (cost) and for parameters and ν were used.

The best results obtained with both -SVR and ν-SVR are shown in Table 5 reveling slightly better results when compared with those obtained with the fitting neural network and better than those obtained with MLP's and the traditional regression algorithms. In all the previous experiments it was tried to estimate G 0 as a function of I D , E D and K D . Following the same strategy as Cruz [START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF], other experiments were performed which results are presented in Table 6. In these experiments G 0 was estimated only as a function of two of the three available variables. As last option E D was used as a single input variable but the results were worst than those obtained using more input variables (due to lack of space they are not shown here). In these experiments only SVR's were used, as they presented the best results in the previous experiments. After the last set of experiments it can be seen that the best results are still those presented in Table 5, the ones that were obtained by using the three available parameters K D , I D and E D .

Type G0 = f (ED, KD) G0 = f (ID, ED) G0 = f (ID,
The final experiments were conducted splitting I D values in coherent groups, namely those related with clay, silt and sandy soils, respectively represented by the following I D intervals: I D < 0.6, 0.6 ≤ I D < 1.8 and I D ≥ 1.8. This operation made a partition of the original 860 data set elements in subsets with 449, 259 and 152 elements respectively. Experiments with these data sets (subsets) were performed using only SVR algorithms since they were the ones with better results in the previous experiments. Results are presented in Table 7. Although these final results can be pointed as very interesting, specially with -SVR, it can not be concluded that these are better than those obtained with the complete data set (Table 5) once there are one larger subset with a result slightly higher and two smaller subsets with a slightly lower result (remember I D -values distribution in Figure 1). This indicates that the methodology can produce interesting results in some particular cases, deserving a more detailed study. Nevertheless from a geotechnical point of view it may be considered a very promising result.

Subsets G0 = f (ID,

Conclusions

A new approach was applied to predict maximum shear modulus by DMT using Neural Networks. Based on performed experiments it is possible to outline the following considerations:

-Neural Networks improve the current state-of-the-art in terms of G 0 prediction through DMT intermediate parameters.

-The results show that, in general, NN lead us to much smaller medians, equivalent means and higher standard deviations in respect to relative errors, when compared to traditional approaches. -Regarding the problem characteristics the SVR approach gives the best results considering the median as the main quality measure as discussed earlier. -The Fitting NN seems to be the most robust in terms of the three quality control parameters used, improving the results of traditional approaches. -The unbalanced data distribution, regarding the I D partition, postpone a final conclusion about the improvement of model quality to the availability of a more balanced sample.

Fig. 1 .

 1 Fig. 1. Sample values for ID,ED,KD and G0

Table 1 .

 1 Sample statistical measures.

		ID	ED	KD	G0
	min	0.0507	0.3644	0.9576	6.4300
	max	8.8143	94.2600	20.5000	529.2000
	median	0.5700	13.4450	3.5750	77.9100
	mean	0.9134	18.8282	4.9161	92.5165
	std	1.0739	18.8264	3.6079	69.6096

Table 2 .

 2 The results obtained with non-linear regression[START_REF] Cruz | Modelling geomechanics of residual soils by DMT tests[END_REF].

		Type	Hidden neurons	G0 = f (ID, ED, KD)
		Quasinewton	50	0.20/0.38(0.72)
	NN	Conj.Grad.	100	0.19/0.30(0.38)
		SCG	40	0.20/0.28(0.33)

Table 3 .

 3 The results (Median/Mean(std)) obtained with Multi Layer Perceptrons.

Table 4 .

 4 The results (Median/Mean(std)) obtained with Bayesian, RBF and Fitting neural networks.

		Type	Hidden neurons	G0 = f (ID, ED, KD)
		MLP-Bayesian	20	0.20/0.29(0.30)
	NN	RBF	200	0.20/0.31(0.39)
		Fitting	60	0.17/0.27(0.29)

Table 5 .

 5 The results (Median/Mean(std)) obtained with Support Vector Regression.

	Type	Cost/ (ν) G0 = f (ID, ED, KD)
	-SVR	200/0.1	0.16/0.27(0.43)
	ν-SVR	200/0.8	0.16/0.27(0.41)

Table 6 .

 6 The results (Median/Mean(std)) obtained with Support Vector Regression using only two input variables.

	KD)

Table 7 .

 7 The results (median/mean(std)) obtained using the subsets of the original data set and for the estimation of G0 as a function of ID, ED and KD.

	ED, KD)


Nick F Ryman-Tubb
email: n.ryman-tubb@surrey.ac.uk

Paul Krause
email: p.krause@surrey.ac.uk

Neural Network Rule Extraction to Detect Credit Card Fraud

Keywords: fraud detection, credit card, neural applications, rule extraction, neuro-symbolic, data mining

Neural networks have represented a serious barrier-to-entry in their application in automated fraud detection due to their black box and often proprietary nature which is overcome here by combining them with symbolic rule extraction. A Sparse Oracle-based Adaptive Rule extraction algorithm is used to produce comprehensible rules from a neural network to aid the detection of credit card fraud. In this paper, a method to improve this extraction algorithm is presented along with results from a large real-world European credit card data set. Through this application it is shown that neural networks can assist in mission-critical areas of business and are an important tool in the transparent detection of fraud.

Introduction

The Oxford dictionary defines fraud as, "wrongful or criminal deception intended to result in financial or personal gain". Fraud can involve the use of stolen credit, debit, fuel or gift cards, opening a bank account through identity theft, concealing the source of illegally or criminally received money, making a wrongful or exaggerated claim for medical expenses through insurance or social security claims, etc. Fraud can be found in on-line shopping, computer gaming, telecommunications, banking, social security, tax evasion, customs fraud, insurance claims, payments and financial services, to the illegal sale of endangered species. The list of crimes and their targets is disappointingly long.

Public perceptions of fraud are often tempered by a belief that it is a "white-collar" crime which targets the wealthy, government and big business and is of less concern as the effects are cushioned for the victim [START_REF]Drawing conclusions about financial fraud: crime, development, and international co-operative strategies in China and the West[END_REF]. However, violent criminals are increasingly moving into fraud [2] so that fraud can involve the threat of violence including murder. In the USA, the fear of fraud supersedes that of terrorism, computer and health viruses and personal safety [START_REF]Research shows economic crisis increases Americans' fears about fraud and ID theft[END_REF]. In the UK the Attorney General described fraud as, "second only to drug trafficking in causing harm to the economy and society." [START_REF] Button | The fraud review and the policing of fraud: Laying the foundations for a centralized fraud police or counter fraud executive?[END_REF]. Today, the proceeds from fraud are paying for organised crime, drug smuggling and terrorism [START_REF] Everett | Credit card fraud funds terrorism[END_REF]. It was estimated that in the UK fraud cost £30bn in 2010, that is 1.3% of the entire UK economic output [6]. In the US, the fraud cost was estimated at $994bn in 2008 [START_REF]Report to the Nation on Occupational Fraud and Abuse[END_REF]. Whatever country, it seems that fraud is large and there is a human cost for each individual act of fraud [START_REF]European Healthcare Fraud & Corruption Network: The Human Cost of Fraud[END_REF].

The detection of fraud is a complex scientific and business challenge. The datasets are large and therefore need to be properly sampled as it becomes computationally impractical to use the natural population. Real world transactional data is noisy, unbalanced, computationally expensive to maintain and highly dimensional. Such data is skewed, has uneven distributions and contains a mixture of symbolic and continuous variables. The dataset is sparse -there is a large quantity of transactional data which contain only a small number of example frauds. Often researchers are unable to report on real fraud datasets due to their sensitivity with many published papers using either synthetic or small datasets; where there is no reason to believe that the conclusions drawn will hold true when they are scaled-up to large, real-world applications [START_REF] Aleskerov | CARDWATCH: a neural network based database mining system for credit card fraud detection[END_REF], [START_REF] Rong-Chang | Personalized Approach Based on SVM and ANN for detecting credit card fraud[END_REF].

Payment Card Fraud

One common type of fraud is payment card fraud -this is the criminal act of deception through the use of a physical plastic card or card information without the knowledge of the cardholder. When a transaction takes place, the details of that transaction are processed by the acquiring bank for authorisation. It is reported that in 2009, the USA total card fraud losses cost banks and merchants $8.6 billion [START_REF] Crosman | Card fraud costs U.S. payment providers $8.6 billion per year[END_REF] and in the UK £609.9 million [START_REF]Card fraud facts and figures[END_REF]. To detect this fraud, organisations use a range of methods; manual methods and some form of automated Fraud Management System (FMS) [START_REF] Cybersource | Seventh annual UK online fraud report: trends, key metrics[END_REF]. The FMS is often a rule-based system that stores and uses knowledge in a transparent way and is easy for a Fraud-Analyst to modify and interpret. Rules provide a convenient mechanism for explaining decisions. However, the generation of comprehensible rules is an expensive and time-consuming task, requiring a high degree of skill both in terms of the developers and the Fraud-Analysts concerned. The performance of the FMS is dependent upon the skill of the Fraud-Analyst and how past data and events are interpreted by them.

An alternative to rules is the use of a learning approach that does not require expert knowledge but learns from examples given in transaction data. A model is formed that can then be used on new transactions to make a decision. The ability of such a model to generalise is fundamental. One such method is to use supervised neural networks that have been widely used to learn from fraud data [START_REF] Aleskerov | CARDWATCH: a neural network based database mining system for credit card fraud detection[END_REF], [START_REF] Brause | Neural data mining for credit card fraud detection[END_REF], [START_REF] Ghosh | Credit card fraud detection with a neural network[END_REF].

Supervised neural networks are essentially a collection of large numbers of realvalued parameters (weights, etc.) with no obvious method to determine their meaning. The knowledge is represented by the distributed weights between the connections, the threshold values and activation functions. The uptake and use of neural networks has been hampered by their "black box" nature and the requirement for often proprietary software to implement the neural network model to be deployed on a server within a mission-critical part of the business -a serious barrier-to-entry in the automated use of FMS.

Mission Critical Application

Once a transaction is made a decision has to be taken to accept or refer/decline the payment. This authorisation process is part of a real-time payments system, which means that the FMS is mission critical -it's failure will cause damage to the business in terms of both money and reputation. It is for this reason that many businesses are reluctant to deploy a FMS that is based on a neural network, where the exact method of fraud detection is kept hidden and they are reliant on their supplier. Many businesses require a transparent approach to fraud detection so that both Fraud-Analysts and management teams can understand exactly why a decision to refer/decline is being made.

A practical solution to promote the widespread use of neural networks within FMS is to use them off-line. One approach is to use the neural network to create evolving models that automatically learn from patterns in transaction data of payment cards that can then be used to extract knowledge in the form of symbolic rules. The rules can be formatted into human-readable ANSI SQL statements that can then be deployed within an existing live Fraud Server environment. This will allow existing FMS systems (that support SQL) to be used and reduce the need for expensive and time consuming creation of manual rules and reliance on Fraud-Analysts.

Rule Generation

There are a number of rule induction algorithms that can be used to generate rules from a dataset. One common method is the decision tree -based on information gain [START_REF] Quinlan | Induction of decision trees[END_REF]. In applications where the dataset is large and contains noise, such as fraud detection, this approach and others have been found to generate a large number of rules where each rule has many conditions and are unfortunately difficult to understand. A key objective for a FMS is to use as few rules as possible -as a "global" view of the fundamental fraud factors is generally preferred over actual accuracy. Neural networks have shown considerable promise in terms of accuracy with generalisation for fraud detection. There are two key approaches to rule extraction from a neural network: (1) Decompositional, (2) Pedagogical. A decompositional approach is where rules are created using heuristic searches that guide the process of rule extraction [START_REF] Towell | Extracting refined rules from knowledge-based neural networks[END_REF], [START_REF] Setiono | Extracting rules from neural networks by pruning and hidden-unit splitting[END_REF] by decomposing the neural network architecture and therefore produce rules that represent the internal architecture. Since there is no reason for individual neurons to represent a recognisable "concept", the extracted rules are often not sufficiently comprehensible. A pedagogical approach is where a set of global rules are extracted from the neural network in terms of the relationships between only the inputs and the outputs [START_REF] Craven | Using sampling and queries to extract rules from trained neural networks[END_REF], [START_REF] Thrun | Advances in Neural Information Processing Systems-Extracting rules from artificial neural networks with distributed representations[END_REF]. These individual rules are then combined into a set of rules that describes the neural network as a whole. The main difficulty with this approach is that the size of the search space is large so that a straightforward search is not practical in real-world problems, where even a small number of input neurons (fields) mean an unrealistic level of computing power is required.

The Sparse Oracle-based Adaptive Rule (SOAR) extraction algorithm is detailed in [START_REF] Ryman-Tubb | SOAR -Sparse Oracle-based Adaptive Rule extraction : knowledge extraction from large-scale datasets to detect credit card fraud[END_REF] and uses sensitivity analysis to avoid the exhaustive decision boundary searches of other neural rule extraction algorithms e.g. [START_REF] Barakat | Learning-based rule-extraction from support vector machines[END_REF], [START_REF] Barakat | Eclectic rule-extraction from support vector machines[END_REF], [START_REF] Engelbrecht | Engineering applications of bio-inspired artificial neural networks-Rule improvement through decision boundary detection using sensitivity analysis[END_REF]. The SOAR algorithm shown in Figure 1, is independent of the neural structure and here uses a standard Multi-Layer Perceptron (MLP) as a Neural "Oracle" which was chosen as a classifier to produce good generalisation and noise-tolerance. Sparse fraud examples are used as the initial search space "seeds" -since (randomly) locating fraud points on the decision boundary in a large search space would be otherwise inefficient. SOAR has the following steps:

1. Continuous valued fields are converted into discreet literals using a well known clustering algorithm called ART2 [START_REF] Carpenter | ART2 : Self-organization of stable category recognition codes for analog input patterns[END_REF] that is applied to group together similar fraud examples. The algorithm creates a new cluster when the input does not belong to the cluster that was determined as most probable, based on a user defined parameter. The Euclidean distance measure is used for measuring similarities between the input pattern and the exemplars in each cluster. 2. The trained MLP neural system is used as an oracle, to interrogate. 3. Prototypes are formed by grouping similar examples from the fraud dataset class to reduce complexity. 4. The prototype generated is "expanded" to cover the largest area on the decision boundary that continues to represent a single class e.g. {genuine, fraud}. This is accomplished using the neural network as an oracle where the binary digits in each literal in the prototype are sequentially activated e.g., {1,0,…0}, {0,1,…0} until {0,0,…1} and the class membership determined by the oracle. In the case of a discretised numeric type each b i represents a contiguous real range and so the search is simplified by only having to activate each in turn, ({b 1 ,b 2 ,b m }, ⊕)=1. Only fields that represent continuous values are expanded since it would make no sense to expand an unordered nominal field. The aim is to capture the generalisation of the neural model while covering as much of the space contiguously as is possible. This is then used to create a propositional rule as a list of antecedents as a combination of the expanded input fields. 5. The rules are optimised to avoid overlap, duplication and redundancy.

This step includes rule pruning to remove rules that have a large estimated error and those which never "fire" (have zero coverage) for the training set.

A false-positive ratio is calculated for each rule using the Training

Dataset. Those rules that have a large ratio (i.e. produce a large number of false-positives compared to the correct number of frauds detected) are removed.

Since the fraud data examples are sparse, the above procedure is computationally efficient. This approach will approximate the classification capability of the neural network by creating a set of rules. The optimisation steps 5 and 6 above aim to improve comprehensibility by reducing the number of rules and approximation errors.

Extended SOAR Extraction

In [START_REF] Ryman-Tubb | SOAR -Sparse Oracle-based Adaptive Rule extraction : knowledge extraction from large-scale datasets to detect credit card fraud[END_REF] the algorithm only extracted rules for the fraud class. The algorithm has been extended here, called SOAR II, so that steps 3 and 4 are repeated to extract rules for a sample of the genuine class. For those points that overlap/intersect within the two sets of rules, false positives will be generated and the fraud rule is "trimmed" so as to remove the overlap. To compare the performance of the two versions of SOAR extraction experiments were completed using the same dataset with both approaches.

Experimental Results

A dataset was provided by a large European company that processes a high volume of transactions per day. The fraud class priors are highly skewed. It has been shown that such a class imbalance may not allow learning algorithms to create good models [START_REF] Chawla | SMOTE: Synthetic Minority Over-sampling TEchnique[END_REF] so sampling has been used. This sampled dataset was created from a natural population of transactions over a fixed time period so as to be representative of the natural population containing all of the known fraud, shown in Table 1.

Pre-Processing and Sampling

The pre-processing was designed to enhance the modelling characteristics of the provided dataset through the following steps:

 All records checked for incomplete or missing data. These records are ignored if found.  Additional feature fields (such as average spend in previous 30-days) are calculated and added to each transaction record.  Attempts are made to identify outliers in the data and to remove them.  Assessing the relevance of each field (to ensure maximum entropy).  Identifying redundancy in the data.  Relevance determination.  Redundancy removal.  Frequency histograms of numeric data are created to ascertain which variables could benefit from transformation.  Frequency analysis of symbolic data is undertaken to ascertain which low frequency variables could benefit from being grouped (into SQL tables).  Re-weighting and splitting.

The pre-processing splits the transaction data into two datasets, shown in Table 2.

Table 2. Sampled Datasets.

Training Dataset Consists of a randomly selected proportion of the transaction data where the class is known. This information is used by the training algorithm to train the neural network to learn the relationships within the data.

Validation Dataset

Consists of a sample of data where the outcome is known and is used only once the model has been completed to validate the accuracy of the model.

Neural Network Training

A three layer MLP architecture was selected as is was the common choice in the previously mentioned research and so is representative of the range of classifiers that could be used. All neuron values were in the range [0,1] using a standard sigmoidal activation function. The training of the neural network is an iterative process that requires considerable processing power to complete where the more transactions presented the longer the algorithm will take. For this reason a Training Dataset was created as a sub-sample shown in Table 3. The Training Dataset used a small subsample of 40,000 transaction records with 2,800 example fraud transactions. The remaining 1,200 fraud transactions were reserved as "unseen" in a Validation Dataset to produce performance metrics given in Table 4. The MLP was then trained using this Training Dataset by interleaving and repeating the fraud examples with the genuine examples so that they were re-balanced. The Conjugate Gradient Descent (CGD) [START_REF] Fletcher | A rapidly convergent descent method for minimization[END_REF] training algorithm was then used. Using a Constructive Algorithm [START_REF] Hirose | Backpropagation algorithm which varies the number of hidden units[END_REF] the number of hidden neurons in the hidden layer was automatically determined by starting at a single neuron and then dynamically growing the layer by adding additional neurons until the stopping criteria is reached; in this case, the final MLP had 60 discreet inputs, 10 hidden and 1 output neuron.

Rule Extraction

Once the neural network was trained, the SOAR and SOAR II extraction algorithms were applied to the Training Dataset and two rule sets created. A confusion matrix was used to evaluate the performance of both the trained neural network and the extracted rule set using the Validation Dataset. The metrics used for calculating accuracy and precision are given in equations ( 1), (2) and (3).

ݕܿܽݎݑܿܿܽ =

ܶܲ + ܶܰ ܶܲ + ܰܨ + ܲܨ + ܶܰ .

(1)

݅ܿ݁ݎ ݅ݏ ݊ ௨ = ܶܲ ܶܲ + ܲܨ . (2) ݅ܿ݁ݎ ݅ݏ ݊ ௨ௗ = ܲܨ ܶܲ + ܲܨ . (3)
While the standard evaluation metrics are important, in the practical application of the approach these need to be projected into the real world. Businesses that deploy a FMS generally have a case management system to review the fraud alerts generated by the rules and have a maximum number of alerts that they can reasonably review in a day. Therefore, performance needs to be projected to the natural population with reference to the number of alerts per day:

ܵܽ݉ ݈ ݅ ݊݃ ܵ݁ݐܴܽ = ݅ݐܿܽݏ݊ܽݎܶ# ݈ܽݎݑݐܽܰݏ݊ ݅ݐܿܽݏ݊ܽݎܶ# ݉ܽܵݏ݊ ݈ ݁ . ( 4
)
ܵ݁݃݉ ݅ݐܽݐ݊݁ ݊ ܴ݁ݐܴܽ ௧ = ݅ݐܿܽݏ݊ܽݎܶ# ݅ܽݎܶݏ݊ ݊݅ ݊݃ ݅ݐܿܽݏ݊ܽݎܶ# ݈ܸܽݏ݊ ݅ ݅ݐܽ݀ ݊ .

(

) 5
݅ݐܿܽݏ݊ܽݎܶ ݊ ܲܨ ݅ݐܴܽ  = ൫FP. ܵ ݃ . ܴ ݐ ൯

TP .

(6)

݅ݐܿܽݏ݊ܽݎܶ ݊ ݈ܽ ݕܽ݀ݎ݁ݏݐݎ݁ ܶ ௗ = ൫.ௌ .ோ ൯ା(.ோ ) ௬௦ ௗ . ( 7
)
݅ݐܿܽݏ݊ܽݎܶ ݊ ݈ݐܿ݁ݎݎܿ ݕ ݈ܽ ݀݁ݐݎ݁ ݕܽ݀ݎ݁ ܶ ௗ = ܶ ௗ ܲܨ . (8)
Here, S g is calculated as 53 (4), R t is 0.04 [START_REF] Everett | Credit card fraud funds terrorism[END_REF] and the resulting performance metrics are given in Table 4, where a "card" is a sequence of transactions linked by a single card number. The number of antecedents in each rule is a good measure on human comprehensibility; shorter and simpler rules are more straightforward to understand than complex rules.

SELECT * FROM [TRANSACTIONS] WHERE RECORDTYPE IN (SELECT A FROM SA) AND COUNTRY IN (SELECT A FROM L1) AND CODE IN (SELECT A FROM CODE1) AND PRDCODE=10 AND NETW IN (SELECT A FROM ID1) AND TRANSTYPE=5 AND IND=1;

Analysis

From Table 4, SOAR II has extracted three fewer rules than the original algorithm with a similar accuracy but at the expense of the fraud precision; those rules that have been omitted in SOAR II were generated from genuine transactions that completely overlap the fraud transactions in the Training Dataset -they are therefore impossible to classify into a single class as they are indistinguishable from each other. The rules in both cases outperform the neural network in terms of accuracy. It is suggested that the rules are able to outperform as the boundaries located are crisp boundaries so that the (inaccurate) generalisation of the neural model is not captured in the extraction. This tends to create rules that do not cover the decision boundary correctly as evidenced by the reduced precision which tends to reduce the false-positive rate suffered by the neural network. The false-positive rate is further reduced by pruning those rules that exhibit a high rate, reducing precision but contributing to improved accuracy. There is no such equivalent process for the neural network. The Fraud-Analysts found that the SOAR II rules were easier to understand and identified both already known patterns as well as previously unknown patterns of fraud.

Conclusions

Accurate rules can be extracted from a trained neural network that can then be inspected and applied in a live FMS environment. This removes the fear of applying neural networks within a mission critical part of a business. The rules are easy to understand and can be produced rapidly, reducing the workload on Fraud-Analysts. SOAR II extraction used on a regular basis has the promise of enabling a more dynamic FMS that can respond to the changing patterns of fraud. It is anticipated that the SOAR II extraction approach can be applied to many other areas of fraud detection where transparency of the decision is important.

Figure 1 .

 1 Figure 1. SOAR Extraction Framework.

  IF RECORDTYPE=(SA) AND COUNTRY=(L1) AND CODE=(CODE1) AND PRDCODE=(10) AND NETW=(ID1) AND TRANSTYPE=(5) AND IND=(1) THEN FRAUDIn SQL:

Table 1 .

 1 European Credit Card Fraud Dataset for Natural Population.

	Feature	Total
	Transactions	60m
	Cards	8m
	Number of Fraud Transactions	4,000
	Value of Fraud Transactions	€1m
	A priory transaction fraud	1:15,000

Table 3 .

 3 Dataset Distribution.

		Total Sample	Training Dataset	Validation Dataset
	Total Transaction Records	1,100,000	40,000	4%	1,060,000 96%
	Containing Fraud Records	4,000	2,800 70%	1,200 30%
	Value of Fraud Records	€1,000,000	€703,000		€297,000

Table 4 .

 4 Evaluation using the Validation Dataset only.For reasons of space, just the top rule from the SOAR II extraction algorithm is presented -this single rule captures over 30% of all frauds in the Validation Dataset. Separate SQL tables of values are created during the pre-processing clustering/grouping stage, e.g. L1 could contain {France, Germany, Spain, Italy, Netherlands}, meaning true if any one of the countries listed is in the transaction field.

	Neural SOAR I	SOAR II

Acknowledgments. This work was supported in part by Retail Decisions Europe Ltd.


Christos M Dimitrakopoulos

Konstantinos A Theofilatos

Efstratios F Georgopoulos

Spyridon D Likothanassis

Athanasios K Tsakalidis

Seferina P Mavroudi
email: mavroudi@ceid.upatras.gr

Adaptive Filtering Techniques Combined with Natural Selection-Based Heuristic Algorithms in the Prediction of Protein-Protein Interactions

Keywords: Protein-Protein Interactions, Adaptive Filtering, Genetic Algorithms. Particle Swarm Optimization, Least Squares Algorithm, Recursive Least Squares, Kalman Filtering

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

Introduction

In each living cell of the human organism, a variety of protein interactions take place. In recent years, researchers have tried to approach the problem of predicting all possible protein interactions in the human organism by implementing different computational techniques. At the beginning, most of them were based on the analysis of a sole feature, indicative of interaction between two proteins. Several examples of such features are features concerning the genomic sequence of the genes-generators of the reference proteins, features concerning the structure of the reference proteins, features concerning the sequences of the references proteins and many others [START_REF] Bock | Predicting protein-protein interactions from primary structure[END_REF][START_REF] Chua | Exploiting indirect neighbors and topological weight to predict protein function from protein-protein interactions[END_REF]. The most recent computational approaches use various features as inputs for their classifiers in order to take advantage of all the available information [START_REF] Chen | Prediction of protein-protein interactions using random decision forest framework[END_REF][START_REF] Fariselli | Prediction of protein-protein interaction sites in heterocomplexes with neural networks[END_REF].

In the present paper, we applied to the problem of PPI prediction several adaptive filtering techniques combined with the most contemporary heuristic methods which are based in the natural selection process. From this combination, a novel computational framework has been formed for the creation of a mathematical equation that gives the best classifier. We consider three classical parameter estimation algorithms (LMS, RLS and Kalman Filter) and one gain adaptation algorithm (IDBD). The gain adaptation algorithms have been shown to perform comparably to the best algorithms (Kalman and RLS), but they have a lower complexity [START_REF] Diniz | Adaptive Filtering: Algorithms and Practical Implementation[END_REF].

Concerning the heuristic methods, each adaptive filtering technique was combined with a genetic algorithm and a Particle Swarm Optimization (PSO) heuristic technique. Genetic algorithms [START_REF] Holland | Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence[END_REF] and PSO [START_REF] Kennedy | Particle Swarm Optimization[END_REF] are the most contemporary heuristic methods and their process is based on the principle of natural selection. They are implemented in order to find the best subset of mathematical terms that constitutes the optimal mathematical model that can be used as the optimal classifier in the adaptive filtering techniques. In the RLS and LMS algorithms, the forgetting and convergence factors are been optimized respectively through the heuristic algorithms used.

In [START_REF] Theofilatos | A new hybrid method for predicting protein interactions using Genetic Algorithms and Extended Kalman Filters[END_REF], genetic algorithms were firstly combined with Extended Kalman Filters for the problem of predicting PPIs. That method has been demonstrated to achieve higher classification performance than the PIPS naive Bayesian method [START_REF] Scott | Probabilistic prediction and ranking of human protein-protein interactions[END_REF] on the same dataset. In the terms of the current paper, the research has been extended by implementing different adaptive filtering techniques and heuristic algorithms in order to establish a novel classification framework. The proposed methodology was tested with a commonly used dataset and the advantages and disadvantages of each method selection are discussed. The mathematical equation that produced the best classifier in terms of classification performance and interpretability was presented and some first biological hypotheses about the classification model were derived.

Materials and Methods

Dataset

The dataset used in this paper was created by retrieving material from the PIPS database [START_REF] Scott | Probabilistic prediction and ranking of human protein-protein interactions[END_REF] which contains information for about 17.643.506 human protein interactions. The positive samples in the dataset were the 16536 known human protein interactions extracted from Human Protein Reference Database (HPRD) [START_REF] Mishra | Human protein reference database 2006 update[END_REF]. The negative dataset was built using randomly chosen protein pairs provided by the PIPS Database. A 1:1 ratio of positive and negative samples was selected in our original dataset, in order to force our classifiers to achieve higher sensitivity and discover a large number of protein-protein interactions. The dataset was split into two equal subsets, training and testing, keeping the 1:1 ratio between the positive and negative samples.

For every protein pair the following features as in [START_REF] Scott | Probabilistic prediction and ranking of human protein-protein interactions[END_REF] were used as inputs:  Expression: Gene expression profiles from 79 physiologically normal tissues (Pearson correlation of co-expression over all conditions).  Orthology: Interactions of homologous protein pairs from yeast, fly, worm and human. The similarity function used is the InParanoid score function [START_REF] O'brien | Inparanoid: a comprehensive database of eukaryotic orthologs[END_REF].  Fully Bayesian combination of the following features:

o Localization: Here PLST predictions are used [START_REF] Scott | Predicting subcellular localization via protein motif co-occurrence[END_REF]. o Domain co-occurrence: Chi-square score of co-occurrence of domain pairs. o Post-translational modifications (PTM) co-occurrence: PTM pair enrichment score has been calculated as the probability of cooccurrence of two specific PTMs in all pairs of interacting protein pairs divided by the probability of occurrence of both of these PTMs separately.  Transitive: This is a module that considers local topology of the underlying network predicted using combinations of above features and it works on the premise that a pair of proteins is more likely to interact if it shares interacting partners. All feature values were normalized in a range from 0 to 1.

Adaptive Filtering Techniques

The adaptive filtering techniques used in this survey are the least mean squares (LMS) [START_REF] Widrow | Adaptive switching circuits[END_REF], the recursive least squares (RLS) [START_REF] Goodwin | Dynamic System Identification: Experiment Design and Data Analysis[END_REF], the Extended Kalman filtering method (EKF) [START_REF] Welch | An Introduction to the Kalman Filter[END_REF] and the incremental delta-bar-delta algorithm (IDBD) [START_REF] Sutton | Adapting bias by gradient descent: An incremental version of delta-bar-delta[END_REF].

The LMS algorithm is by far the most widely used algorithm in adaptive filtering because of its low computational complexity, proof of convergence in stationary environment and stable behavior when implemented with finite-precision arithmetic.

To guarantee the convergence of LMS it is good practice to set the learning rate in the range,

, where λ max is the largest eigenvalue of R = E[x(k)x T (k)], where x(k) is the vector of the input variables.

Least-squares algorithms aim at the minimization of the sum of the squares of the difference between the desired signal and the model filter output. When new samples of the incoming signals are received at every iteration, the solution for the least-squares problem can be computed in recursive form resulting in the recursive least-squares (RLS) algorithms.

Kalman filter is a set of mathematical equations which constitutes an efficient means to estimate the state of a process by minimizing the mean of the square error. The Extended Kalman Filter (EKF) is the nonlinear version of the Kalman Filter, and its goal is to approach the situation where the process to be estimated or the measurement relationship to the process is nonlinear. The specific equations for the time and measurement updates are divided into two groups. The time update equations and the measurement update equations.

Gain adaptation algorithms [START_REF] Sutton | Adapting bias by gradient descent: An incremental version of delta-bar-delta[END_REF] implement a sort of meta-learning in that the learning rate is adapted based on the current inputs and on the trace of previous modifications [START_REF] Ratnaweera | Prediction Update Algorithms for XCSF: RLS, Kalman Filter, and Gain Adaptation 18[END_REF]. The incremental delta-bar-delta (IDBD) uses a different adaptive learning rate for each input. This action can lead to the improvement of the system, when some of the inputs are irrelevant. In IDBD, each element of the k i (t) of the gain vector k(t) is computed separately.

Heuristic Methods

Heuristic methods are implemented to fasten the process of finding a "good enough" solution to a problem, where the usage of an exhaustive search is impractical. Because of the large search space of our problem (2 113 possible solutions), 2 classical heuristic methods were used, genetic algorithms and particle swarm optimization.

Genetic Algorithms (GAs) [START_REF] Holland | Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence[END_REF] are search algorithms inspired by the principle of natural selection. They are useful and efficient when the search space is big and complicated or when there is not any available mathematical analysis of the problem. A population of candidate solutions, called chromosomes, is optimized through a number of evolutionary cycles and genetic operations, such as crossovers or mutations.

In our approach, a simple GA was used where each chromosome comprises term genes that encode the best term subset and parameter genes that encode the best choice of parameters for each adaptive method used. The term genes are binary genes with the value of 1 indicating that this specific term should be included in the final classification model and the value of 0 indicating the opposite. The parameters which are optimized using GA are the convergence factor μ for the LMS adaptive filtering method and the forgetting factor λ for the RLS adaptive filtering method.

For the genetic algorithm used in our hybrid methodology, the one-point crossover and the mutation operators were used. Crossover and mutation probabilities for the GA where set to 0.9 and 0.01 respectively. Crossover is used in hope that new chromosomes will have good parts of old chromosomes and maybe the new chromosomes will be better. However it is good to leave some part of population survive to next generation. This is the reason a high (but not equal to one) crossover probability was used. As already mentioned, mutation is made to prevent falling GA into local extreme, but it should not occur very often, because then GA will in fact change to random search. That is the main reason why a small mutation probability was applied. Furthermore, roulette selection was used for the selection step and elitism to raise the evolutionary pressure in better solutions and to accelerate the evolution.

The size of the initial population was set to 30 chromosomes after experimentation in the training dataset. The termination criterion is the maximum number of 100 generations to be reached combined with a termination method that stops the evolution when the population is deemed as converged. The population is deemed as converged when the average fitness across the current population is less than 5% away from the best fitness of the current population. Specifically, when the average fitness across the current population is less than 5% away from the best fitness of the population, the diversity of the population is very low and evolving it for more generations is unlikely to produce different and better individuals than the existing ones or the ones already examined by the algorithm in previous generations.

The particle swarm optimization (PSO) [START_REF] Kennedy | Particle Swarm Optimization[END_REF] algorithm is another population based heuristic search algorithm based on the simulation of the social behavior of birds within a flock.

In our approach PSO searches the best mathematical term's subset for building the optimal classifier and the best parameters for each adaptive filtering method. Thus its particle is consisted of 113 variables indicating whether to include a term in the subset or not, plus one variable for the parameter selection of the adaptive filtering method. The values of the variables for the term selection range from 0 to 2. Values bigger than 1 indicate that this specific term should be included in the final classification model and the values less than 1 indicate the opposite.

In order to make a fair comparison with the GA search method, the initial population size of the swarm was set to 30 particles and the total number of iterations was set to 100. The termination criteria of convergence was not applied here because of the adaptive search behavior used in our PSO implementation which needs the total number of iterations to be used in order to perform in its full strength and effectiveness.

As mentioned above, the convergence factor μ for the LMS adaptive filtering method and the forgetting factor λ for the RLS adaptive filtering method are been optimized through the GA and the PSO algorithms. According to the theory, the convergence factor μ for LMS algorithm has been selected in the range [0, 1/λ max ], in order to guarantee convergence of the mean square error, where λ max is the largest eigenvalue of the autocorrelation matrix of the input x(k). Concerning RLS algorithm, the forgetting factor λ is chosen such that 2μ = 1 -λ, where μ is the corresponding convergence factor of the LMS algorithm. Hence, the values of the parameter λ are close to 1 and the values of the parameter μ are close to 0.

Evaluation in the Hybrid Methodology

The main idea of our proposed classification method is to use a heuristic method to find a "good enough" subset of terms in order to build the mathematical model for our predictor and then apply an adaptive filtering technique to find its optimal parameters. The search space of our problem consists of 2 113 possible solutions and hence is extremely large. GAs and PSOs are heuristic methods that can deal with large search spaces and do not get trapped in local optimal solutions. The 113 mathematical terms that were used in our method were taken from 3 known nonlinear mathematical models, which are the Volterra Series, the exponential and the polynomial model, as described in [START_REF] Theofilatos | A new hybrid method for predicting protein interactions using Genetic Algorithms and Extended Kalman Filters[END_REF].

In our approach, a simple Genetic Algorithm and PSO were used as heuristic methods. Each chromosome in these methods comprises of genes that encode the best subset of mathematical terms to be used in our classifier. Roulette selection was used for the selection step including elitism to accelerate the evolution of the population in the Genetic algorithm.

The evaluation process in both heuristic techniques is described in the following steps:

 For every individual of the population, use the Adaptive Filtering method (using the training dataset) to compute the best parameters for the subset of terms that arises from the individual's genes.  Use the mathematical model found at the previous step to compute the fitness of the classifier in the validation set.  For every individual compute the output of the following function:

F = a* G m + b*sensitivity + c*specificity -d*terms (2)
, where G m is the geometric mean of sensitivity and specificity. Hence, the proposed method tries to produce one mathematical model that optimizes the combination of sensitivity and specificity measures, sensitivity and specificity alone and simultaneously minimize the number of mathematical terms to be used in the mathematical model.  Rank the individuals from the one having the lower fitness value to the one having the bigger fitness value. Using this ranking in the evaluation procedure, better scaled fitness scores can be achieved.  Give as fitness value for every individual their ranking number. For example, the worst individual will take fitness equal to 1, the next fitness equal to 2 etc.

During the training we used 5-fold cross validation.

Experimental Results

In order to evaluate the performance of each combination of algorithms we applied it to the dataset described in Section 2.1. In order to make a fair comparison between all combinations, their classification thresholds were optimized. For all methods we experimented running them 100 times and in Table 1, the average values for each combination of algorithms are presented. As it is clear from Table 1, Kalman Filtering algorithm when combined with the PSO heuristic method achieves the highest geometric mean, which is the combination of sensitivity and specificity. Hence, we can assume that Kalman with PSO achieves the highest classification performance. However, at the same time it is observed that the optimal mathematical models generated by Kalman with GA are the simplest ones concerning complexity. Particularly, classifiers generated from Kalman with GA algorithm have an average of 12.3 mathematical terms. While IDBD and RLS succeed to achieve higher specificity than Kalman algorithm, they cannot achieve a higher overall classification performance through the geometric mean measure, because of their low sensitivity values.

In general, GAs lead to less complexity classifiers than the classifiers produced by PSO for all possible adaptive filtering algorithms. According to the classification performance, adaptive filtering techniques can be ranked from the best one to the worse as Kalman, RLS, LMS and IDBD. This fact can also be observed at figure 1, where the performance of each algorithm has been depicted based on the evaluation measure, as described in equation [START_REF] Chua | Exploiting indirect neighbors and topological weight to predict protein function from protein-protein interactions[END_REF]. The evaluation measure is a multi-objective measure, as it considers more than one metrics simultaneously. Kalman and RLS perform almost the same, but Kalman achieves a slightly higher accuracy. LMS algorithm is known for its computational simplicity, but RLS algorithm and Kalman filter have higher performance as they focus on the minimization of the mean square error between the desired signal and the model filter output. IDBD, which is a gain adaptation algorithm, fails to compete with the other algorithms. Maybe this happens because the learning rate of the IDBD algorithm is not optimized through the heuristic algorithms, but within the algorithm itself (meta-learning parameter). When dealing with the PPI prediction problem our first objectives are classification performance and interpretability and not computational complexity because of the off-line nature of the problem. Thus, among the local search algorithms explored in the present paper, Extended Kalman filters are the best solution.

Another observation is that the results of the adaptive filtering techniques combined with PSO have larger variance than the corresponding algorithms combined with GA. This means that PSO is more unstable than GA algorithm, and hence, GA has a smoother convergence. Moreover, we can observe that RLS, LMS and IDBD have -in general-larger variance than the Kalman algorithm which is more stable. This fact can be attributed to the parameters that RLS, LMS and IDBD try to optimize through the heuristic algorithm. As a result, they have a larger search space to explore and the complexity of the heuristic problem rises making them difficult to converge using the same iterations. In the Kalman case, both GA and PSO have low variance in the measures' results and PSO achieves to overcome the classification accuracy of the GA. Next, we present the classifier with the higher evaluation value (equation 2), that resulted from the Kalman with GA combination. We observe that there are a lot of nonlinear mathematical terms in the model, like exponential or polynomial, revealing the complexity of the real model and indicating that this classification model cannot be derived with simple classification algorithms. The model attributes higher importance to terms including the third feature x 3 , which is a combination of three distinct features. Furthermore, the third feature dominates the model, as 4 terms out of 9 consist only of the third feature.

(3)

Conclusions

We have studied the influence of different combinations between adaptive filtering techniques and heuristic algorithms on PPI prediction problem in order to produce the optimal hybrid method for this problem. A search space of possible mathematical models is been searched through the heuristic algorithms (GA, PSO) and the optimal parameters for each specific model are been calculated by using an adaptive filtering technique (Kalman, RLS, LMS, IDBD). The mathematical model resulted from the above process constitutes our optimal classifier. All methods were tested in a human protein interaction dataset extracted from the PIPS database.

Our research concluded that the combination of an Extended Kalman Filter combined with PSO -as described in Section 2.3 -gives us the highest classification performance. Moreover, the Extended Kalman Filter when combined with GA produces the less complicated classifiers. In general, PSO heuristic has a larger variance than GA technique, except for the Kalman algorithm where both are stable. In the Kalman algorithm, PSO achieves to overcome the classification accuracy of GA.

In this paper, we presented a new classification framework for predicting PPIs combining heuristic algorithms with adaptive filtering techniques. According to the need of the researchers for high classification accuracy or low complex classifiers, the most convenient combination of algorithms can be used. Finally, we intent to use our method for the prediction of novel human PPIs.

Fig. 1 .

 1 Fig. 1. Evaluation value is the average of the three main measures (geometric mean, sensitivity, specificity), for each combination of algorithms. 1: Kalman GA, 2: Kalman PSO, 3: RLS GA, 4: RLS PSO, 5: LMS GA, 6: LMS PSO, 7: IDBD GA, 8: IDBD PSO.

Table 1 .

 1 Classification performance of all possible combinations between adaptive filtering and heuristic algorithms. Gm is the geometric mean of sensitivity and specificity.

	Algorithm	Gm	Sensitivity Specificity # of terms
	Combination
	Kalman with GA	0.7911	0.7540	0.8302	12.3
	Kalman with PSO	0.7927	0.7568	0.8304	39.3
	RLS with GA	0.7826	0.7291	0.8400	14.6
	RLS with PSO	0.7797	0.7308	0.8319	26.1
	LMS with GA	0.7012	0.6388	0.7698	15.6
	LMS with PSO	0.6845	0.6418	0.7300	38.6
	IDBD with GA	0.6704	0.4976	0.9033	26.6
	IDBD with PSO	0.6694	0.5516	0.8124	37.3

Acknowledgments. This research has been co-financed by the European Union (European Social Fund -ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) -Research Funding Program: Heracleitus II. Investing in knowledge society through the European Social Fund.