Monique Thonnat

Semantic Activity Recognition

Extracting automatically the semantics from visual data is a real challenge. We describe in this paper how recent work in cognitive vision leads to significative results in activity recognition for visualsurveillance and video monitoring. In particular we present work performed in the domain of video understanding in our PULSAR team at INRIA in Sophia Antipolis. Our main objective is to analyse in real-time video streams captured by static video cameras and to recognize their semantic content. We present a cognitive vision approach mixing 4D computer vision techniques and activity recognition based on a priori knowledge. Applications in visualsurveillance and healthcare monitoring are shown. We conclude by current issues in cognitive vision for activity recognition.

INTRODUCTION

This paper is focused on activity recognition. Activity recognition is a hot topic in the academic field not only due to scientific motivations but also due to strong demands coming from the industry and the society; in particular for videosurveillance and healthcare. In fact, there is an increasing need to automate the recognition of activities observed by visual sensors (usually CCD cameras, omni directional cameras, infrared cameras). More precisely we are interested in the real-time semantic interpretation of dynamic scenes observed by video cameras. We thus study spatio-temporal activities performed by mobile objects (e.g. human beings, animals or vehicles) interacting with the physical world.

What does it mean to understand a video ? Is it just to perform statistics on the appearance of images and to recognize an image from a set of already seen images? If we really want to understand the activities performed by the physical objects 2D analysis is not sufficient. We need to locate the physical objects in the 3D real world. The dynamics of the physical objects is a major cue for activity recognition. The computer vision community is very active in the domain of motion detection, mobile object tracking and more recently trajectory analysis. Very often these analyses are performed in the image plane and are thus dependant of the sensor parameters as its field of view, position and orientation. However for reliable activity recognition the dynamics of the physical objects must be computed in the 4D space.

Is there a unique objective interpretation of a dynamic scene? For instance the scenes shown in figures 1 and 2 can be interpreted more or less precisely in function of the a priori knowledge of the observer. In the first case (shown in figure 1) without information on the location of the scene one can recognize an indoor scene where two men are walking together towards a door; a videosurveillance expert knowing the location (a bank agency), its spatial configuration as well as security rules will interpret the same scene as a bank attack 1 INRIA, France, email: Monique.Thonnat@sophia.inria.fr with the unautorized person accessing together with an employee to a fordidden area. In the second case (shown in figure 2) without information on the location of the scene one can recognize a woman standing alone; a medical expert knowing the patient will interpret the same scene as an active elderly preparing a meal in her kitchen. In fact, the interpretation of a video sequence is not unique but it depends on the a priori knowledge of the observer and on his/her goal.

4D APPROACH

We present a cognitive vision approach mixing 4D computer vision techniques and activity recognition based on a priori knowledge. The major issue in semantic interpretation of dynamic scenes is the gap between the subjective interpretation of data and the objective measures provided by the sensors. Our approach to address this problem is to keep a clear boundary between the application dependent subjective interpretations and the objective analysis of the videos. We thus define a set of objective measures which can be extracted in real-time from the videos, we propose formal models to enable users to express their activities of interest and we build matching techniques to bridge the gap between the objective measures and the activity models.

Figure 3 shows the global structure of a videosurveillance system built with this approach. First, a motion detection step followed by a frame to frame tracking is made for each video camera. Then the tracked mobile objects coming from different video cameras with overlapping fields of view are fused into a unique 4D representation for the whole scene. Depending on the chosen application, a combination of one or more of the available trackers (individuals, groups and crowd tracker) is used. Then scenario recognition is performed by a combination of one or more of the available recognition algorithms (automaton based, Bayesian-network based, AND/OR tree based and temporal constraints based). Finally the system generates the alerts corresponding to the predefined recognized scenarios.

For robust semantic interpretation of mobile object behaviour it is mandatory to rely on correct physical object type classification. It can be based on simple 3D models like parallelepipeds [START_REF] Zúniga | Fast and reliable object classification in video based on a 3d generic model[END_REF] or complex 3D human body configurations with posture models as in [START_REF] Boulay | Applying 3d human model in a posture recognition system[END_REF]. Figure 4 shows examples of such postures.

3D MAP

We use 3D maps as a means to model the a priori knowledge of the physical environment captured by the sensors. More precisely the 3D maps contain the a priori knowledge of the empty scenes:

• Video Cameras: 3D position of the sensors, calibration matrix, fields of view,... • 3D Geometry: the geometry of the static structure of the empty scene (for instance the buildings and road structure for outdoor scenes or the walls and doors for indoor scenes) as well as the main static 3D objects (for instance the furniture in indoor scenes) and the 2D zones of interest. This geometry is defined in terms of 3D position, shape and volume. • Semantic information: for each part of the map semantic information is added as its type (e.g. 3D object, 2D zone), its characterics (e.g. yellow, fragile) or its function (e.g. entrance zone, seat).

We can see on figure 5 a 2D map of an indoor flat and on figure 10 two partial views of the 3D map built for monitoring elderly at home. In this map in addition to the main structure of the rooms (walls, doors, etc.), the equipment and the furniture are defined as well as the information related to the sensors.

ACTIVITY MODELLING

In order to express the semantics of the activities a modelling effort is needed. The models correspond to the modeling of all the knowledge needed by the system to recognize video events occurring in the scene. To allow security operators to easily define and modify their models, the description of the knowledge is declarative and intuitive (in natural terms). We propose a video event ontology to share common concepts in video understanding and to decrease the effort of knowledge modelling.

The Video Event Ontology

The event ontology is a set of concepts for describing physical objects, events and relations between concepts:

The physical objects are all the concepts to describe objects of the real world in the scene observed by the sensors. The attributes of a physical object are pertinent for the recognition. These attributes characterize the physical object. There are two types of physical objects: contextual objects (which are usually static and whenever in motion, its movement can be predicted using contextual information) and mobile objects (which can be perceived as moving in the scene and as initiating their motions, without the possibility to predict their movement).

The events are all the concepts to describe mobile object evolutions and interactions in a scene. Different terms are used to describe these concepts and categorized into two categories: state (including primitive/composite state) and event (including primitive/composite event, single/multi-agent event).

A primitive state is a spatio-temporal property valid at a given instant or stable on a time interval which is directly inferred from audiovisual attributes of physical objects computed by low level signal processing algorithms.

A composite state is a combination of states. A primitive event is a change of states. A composite event is a combination of states and events. A single-agent event is an event involving a single mobile object. A multi-agent event is a composite event involving several (at least two) mobile objects with different motions.

Currently this ontology contains 151 concepts used for different applications in video understanding. This ontology is implemented in Protege to be independant of a particular activity recognition formalism.

Activity Models

A formalism for expressing an activity is directly based on the concepts of the video event ontology. A composite event model is composed of five parts: "physical objects" involved in the event (e.g. person, equipment, zones of interest), "components" corresponding to the sub-events composing the event, "forbidden components" corresponding to the events which should not occur during the main event, "constraints" are conditions between the physical objects and/or the components (including symbolic, logical, spatial and temporal constraints including Allen interval algebra operators, and "alarms" describing the actions to be taken when the event is recognized.

Primitive states, composite states and primitive events can be described using the same formalism. Please see [START_REF] Vu | Automatic video interpretation: A recognition algorithm for temporal scenarios based on pre-compiled scenario models[END_REF] and [START_REF] Vu | Automatic video interpretation: A novel algorithm for temporal scenario recognition[END_REF] for more details of the formalism.

ACTIVITY RECOGNITION

The algorithm proposed in [START_REF] Vu | Automatic video interpretation: A novel algorithm for temporal scenario recognition[END_REF] and in [START_REF] Vu | Automatic video interpretation: A recognition algorithm for temporal scenarios based on pre-compiled scenario models[END_REF] enables to process efficiently (i.e. in realtime) a data flow and to recognize pre-defined activities. Alternative approaches based on probabilistic methods [START_REF] Medioni | Activity Analysis in Video[END_REF] or [START_REF] Moenne-Loccoz | Recurrent bayesian network for the recognation of human behaviors from video[END_REF] can also be used. In the following we concentrate on the first approach because it is directly based on the formalism and the ontology presented in the previous section. The video event recognition algorithm recognizes which events are occurring using the primitive video events. To recognize an event composed of sub-events, given the event model, the recognition algorithm selects a set of physical objects matching the remaining physical object variables of the event model. The algorithm then looks back in the past for any previously recognized state/event that matches the first component of the event model. If these two recognized components verify the event model constraints (e.g. temporal constraints), the event is said to be recognized. In order to facilitate complex event recognition, after each event recognition, event templates are generated for all composite events, the last component of which corresponds to this recognized event. For more details see [START_REF] Vu | Automatic video interpretation: A novel algorithm for temporal scenario recognition[END_REF].

APPLICATIONS

This approach has been applied to a large set of applications in visualsurveillance.

Visualsurveillance

A typical example of complex activities in which we are interested is aircraft monitoring (see figure 7 in apron areas . In this example the duration of the servicing activities8 around the aircraft is about one hour and the activities involve interactions between several ground vehicles and human operators.

The goal is to recognize these activities through formal activity models as shown in figure 9 and data captured by a network of video cameras (such as the ones shown in figure 7). For more details, refer to [START_REF] Fusier | Video understanding for complex activity recognition[END_REF]

Healtcare monitoring

In this application the objective is to monitor elderly at home (see figure 10). In collaboration with gerontologists, we have modeled several primitive states, primitive events and composite events. First we are interesting in modelling events characteristic of critical situations such as falling down. Second, these events aim at detecting abnormal changes of behavior patterns such as depression. Given these objectives we have selected the activities that can be detected using video cameras [START_REF] Zouba | Monitoring activities of daily living (adls) of elderly based on 3d key human postures[END_REF]. We have modeled thirty four video events. In particular, we have defined fourteen primitives states, four of them are related to the location of the person in the scene (e.g. inside kitchen, inside livingroom) and the ten remaining are related to the proposed 3D key human postures. We have defined also four primitive events related to the combination of these primitive states: "standing up" which represents a change state from sitting or slumping to standing, "sitting down" which represents a change state from standing, or bending to sitting on a chair, "sitting up" represents a change state from lying to sitting on the floor, and "lying down" which represents a change state from standing or sitting on the floor to lying. We have defined also six primitive events such as: stay in kitchen, stay in livingroom. These primitive states and events are used to define more composite events. For this study, we have modeled ten composite events. In this paper, we present just two of them: "feeling faint" and "falling down".

The model of the "feeling faint" event is shown in figure 4. The "feeling faint" model involves one physical object (one person), and it contains three 3D human posture components and constraints between these components.

CONCLUSION

We have shown a 4D semantic approach for activity recognition of dynamic scene. There are still a lot of open issues among which a full theory of visual data interpretation, reliable techniques for 4D analysis able to deal with changing observation conditions and scene content. From an activity recognition point of view the three main points are the development of shared operational ontologies, of formalisms for activity modelling with good properties such as scalability and learning techniques for model refinement. In particular a large set of learning issues are rised by this 4D semantic approach for instance: learning contextual variations for physical object detection and image segmentation [START_REF] Martin | Learning contextual variations for video segmentation[END_REF], learning the structure of the activity models [START_REF] Toshev | An a priori-based method for frequent composite event discovery in videos[END_REF] or learning the visual concept detectors [START_REF] Maillot | Ontology based complex object recognition[END_REF].

Figure 1 .

 1 Figure 1. A scene with different valid interpretations: two people walking together towards a door or a bank attack with an access to a forbidden area by an unauthorized person and an employee.

Figure 2 .

 2 Figure 2. A scene with different valid interpretations: a person standing in a room or an active elderly preparing a meal in a kitchen.

Figure 4 .

 4 Figure 4. Different 3D models of human body postures

Figure 5 .

 5 Figure 5. Top view of the flat

Figure 6 .

 6 Figure 6. 3D map: the kitchen area and the top view of a flat for monitoring elderly at home

Figure 7 .

 7 Figure 7. Different views of an apron area captured by video cameras for aircraft monitoring

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Activity recognition problem in airport: the main servicing operations around an aircraft (refuelling, baggage loading, power supply, etc...) and the location of the 8 video cameras (in blue)

  CompositeEvent (PersonFeelingFaint, PhysicalObjects( (p: Person) ) Components ( (pStand: PrimitiveState Standing(p)) (pBend: PrimitiveState Bending(p)) (pSit: PrimitiveState Sitting Outstretched Legs(p))) Constraints ((Sequence pStand; pBend; pSit) (pSit's Duration >= 10)) Alarm( AText("Person is Feeling Faint") AType("URGENT")) ) "Feeling faint" model. We have also modelled the "falling down" event. There are different ways for describing a person falling down. Thus, we have modelled the event "falling down" with three models: Falling down 1: A change state from standing, sitting on the floor (with flexed or outstretched legs) and lying (with flexed or outstretched legs). Falling down 2: A change state from standing, and lying (with flexed or outstretched legs). Falling down 3: A change state from standing, bending and lying (with flexed or outstretched legs). An example of the definition of the model "falling down 1" is shown below. CompositeEvent(PersonFallingDown1, PhysicalObjects( (p: Person) ) Components ( (pStand: PrimitiveState Standing(p)) (pSit: PrimitiveState Sitting Flexed Legs(p)) (pLay: PrimitiveState Lying Outstretched Legs(p))) Constraints ( (pSit before meet p lay) (pLay's Duration >= 50)) Alarm (AText("Person is Falling Down") AType("VERYURGENT")) )"Falling down 1" model.

Figure 11

 11 Figure 11 and figure12show respectively the camera view and the 3D visualization of the recognition of the "feeling faint" event.

Figure 11 .Figure 12 .

 1112 Figure 11. Recognition of the "feeling faint" event

Figure 13 .

 13 Figure 13. Recognition of the "falling down" event

  From sensor data to high level interpretation; global structure of an activity monitoring system built with VSIP[START_REF] Avanzi | Design and assessment of an intelligent activity monitoring platform[END_REF].

	images from camera 1	motion detection	tracking frame to frame mobile objects from camera 1	fused tracked mobile objects for the whole scene tracked mobile objects from camera 1	long term individual tracking	AND/OR tree-based scenario recognition automaton-based
	images from camera 2	motion detection	tracking frame to frame	cameras with overlapped FOV fusion	tracking long term group	temporal-constraints based scenario recognition	alerts
	images from camera N	... motion detection	... tracking frame to frame	tracked mobile objects from camera N	long term crowd tracking physical objects	scenario recognition Bayesian-network based scenario recognition
	Figure 3.


come

Conception d'Environnements Informatiques Communautaires pour MAETIC

Céline Quénu-Joiron


Rémi Rhodes

Vincent Vargas

J C Domenge

FORECASTING VOLATILITY IN THE PRESENCE OF LEVERAGE EFFECT

We define a simple and tractable method for adding the Leverage effect in general volatility predictions. As an application, we compare volatility predictions with and without Leverage on the SP500 Index during the period 2002-2010.

Introduction

In this note, we adress the following problem: starting from a quite general symmetrical model of returns, how does one perturb it naturally in order to get a model with Leverage Effect? In formulas (2.1), (2.2) below, we propose a simple framework in this direction. The main motivation for this work is the problem of forecasting volatility , which has applications in many fields of finance: risk management, option pricing, etc... The most interesting feature of the construction is precisely that one can derive from the symmetrical model simple prediction formulas (cf. formula (5.1) below) for the perturbed model.

The paper is organized as follows: Section 2 defines the model in the discrete case and states the main properties of the model: correlation functions, existence of moments,etc... Section 3 defines the model in the continuous case. Section 4 is devoted to the simulation problem. Section 5 adresses the issue of forecasting volatility. Finally, we gather in the appendix the proofs of Section 2,3,4.

The discrete case

2.1. The symmetrical model. We consider the general stationary model for the (log) returns given by: r i = σ i ǫ i , where ǫ = (ǫ i ) i∈Z is a sequence of centered i.i.d. variables such that E[ǫ 2 i ] = 1 and the volatility (σ i ) i∈Z is a random process that we will write under the form:

σ i = σ(γ + X i ),
where X = (X i ) i∈Z is a centered stationary process independent of ǫ with E[X 2 i ] < 1 and γ 2 = 1 -E[X 2 i ]. With these conventions, we get that E[r 2 i ] = σ 2 . In financial applications, X will have long range correlations:

1 Université Paris-Dauphine, Ceremade, Paris, France, rhodes@ceremade.dauphine.fr 2 Université Paris-Dauphine, CNRS, Ceremade, Paris, France, vargas@ceremade.dauphine.fr 3 Capital Fund management, 6, blvd Haussmann 75009 Paris, URL: http://www.cfm.fr/

1 • Power law: E[X i X i+j ] ≈ A/(j + 1) µ . Typically, we will take µ ∈ [0, 1]. • Multifractal (µ → 0): E[X i X i+j ] ≈ λ 2 ln + (T /(j + 1)
). Typically, we will take λ 2 ∈ [0.02, 0.05] and T ≈ 2000.

2.2.

The model with leverage effect. We consider the sequences ǫ, X of the previous section. We introduce two extra parameters β, α respectively standing for the magnitude of the leverage effect and for the inverse of the relaxation time of the leverage effect. We want to consider the following model:

r i = σ i ǫ i , (2.1)
where the volatility (σ i ) i∈Z is a random process that satisfies the following recursive equation:

σ i = σ(γ + X i -β i-1 k=-∞ e -α(i-k) r k ), (2.2)
where

γ 2 = 1 -E[X 2 i ] -σ 2 β 2 e 2α -1
. We introduce the filtration:

F i = σ{(X j ) j i , (ǫ j ) j i-1 }
This leads to the following natural definition:

Definition 2.1. We say that a sequence

(σ i , r i ) i∈Z solution of (2.1), (2.2) is non anticipative if (σ i ) i∈Z is F i -adapted.
We can now state the following existence theorem:

Theorem 2.2. There exists a unique stationary non anticipative square-integrable solution (σ i , r i ) i∈Z of (2.1) and (2.2) if and only if σ 2 β 2 e 2α -1 < 1. In the sequel of this section, we will therefore consider the unique stationary and non anticipative solution of (2.1), (2.2). It is straightforward to compute the following quantities (see appendix):

• Average vol and Variance:

E[σ i ] = γσ and E[r 2 i ] = E[σ 2 i ] = σ 2 .
• Volatility fluctuations:

E[σ i σ i+j ] -E[σ i ] 2 = σ 2 (E[X i X i+j ] + β 2 e 2α -1 e -αj ).
• Leverage correlations:

E[r i σ i+j ] E[r 2 i ]E[σ i+j ] = -β γ e -αj Remark: Assuming the renormalisation relation γ 2 = 1 -E[X 2 i ] -σ 2 β 2 e 2α -1
is not a restriction. Indeed, if we are given a set of coefficients (σ, γ, β, α) and a stationary process X that do not satisfy this relation, it is possible to renormalize the coefficients without changing the solution of (2.1) and (2.2). It suffices to define

c = γ 2 + E(X 2 i ) 1 -β 2 σ 2 e 2α -1 and γ = γ c , σ = cσ, β = β c , X i = X i c .
Then the new set of coefficients ( σ, γ, β, α) and the stationary process X satisfy the renormalisation relation. Moreover, the solution of (2.1) and (2.2) is clearly the same for the datas (σ, γ, β, α, X) or ( σ, γ, β, α, X).

The continuous case

A natural extension of our discrete model to the continuous case is the following. Given a standard Brownian motion B = (B t ; t ∈ R) and a centered square-integrable stationary process X = (X t ; t ∈ R) independent from B, we are interested in finding a stationary solution (σ t ) t∈R to the equations:

dR t = σ t dB t , (3.1)
σ t = σ γ + X t -β t -∞ e -α(t-r) σ r dB r , (3.2)
where α, β, γ, σ are positive parameters satisfying the relation

γ 2 = 1-E[X 2 t ]-σ 2 β 2 e 2α -1 .
Once again, we introduce the filtration:

F t = σ{(X s ) s t , (ǫ s ) s t },
leading to the following notion of non-anticipativity Definition 3.1. We say that a sequence (σ t , r t ) t∈R solution of (3.1)

, (3.2) is non anticipative if (σ r ) t∈R is F t -progressively measurable.
We can now state the following existence theorem: Theorem 3.2. There exists a unique stationary non anticipative square-integrable solution (σ i , r i ) i∈Z of (2.1) and (2.2) if and only if σ 2 β 2 2α < 1. The unique stationary and non anticipative solution of (3.1), (3.2) satisfies:

• Average vol and Variance:

E[σ t ] = γσ, E[(R t -R 0 ) 2 ] = tσ 2 , and E[σ 2 t ] = σ 2 .
• Volatility fluctuations: for s > 0,

E[σ t σ t+s ] -E[σ t ] 2 = σ 2 (E[X t X t+s ] + β 2 2α e -αs ).
• Leverage correlations: for s > h > 0,

E[(R t+h -R t )σ t+s ] E[(R t+h -R t ) 2 ]E[σ t+s ] = - β γ e
-αs e αh -1 h .

Simulation

Since the solution of (2.1) and (2.2) is obtained by the Picard fixed point theorem over a space of stationary processes and since the evolution at time i of σ depends on the whole trajectory of σ between -∞ and i -1, the simulation of such a process is not straightforward. We suggest the following method based on finite dependence. Our purpose is to simulate the process (σ i , r i ) i∈Z between dates 0 and p > 0.

We choose N > 1 and define the mapping T N which maps a sequence u = (u i ) i∈Z to the sequence u = ((T N u) i ) i∈Z by:

(T N u) i = σ(γ + X i -β i-1 k=i-N e -α(i-k) u i ǫ i ).
We define recursively the sequence (σ N,n ) n∈N of random sequences by σ N,0 = (0) i∈Z and ∀n 1, σ N,n = T N (σ N,n-1 ).

We get the following approximation result: Theorem 4.1. Given N 1, the approximation error between (σ N,n i ) i∈Z and the solution (σ i ) i∈Z of equations (2.1) and (2.2) can be estimated by

sup i∈Z E[|σ i -σ N,n i | 2 ] σ 2 e -αN 1 -C 1/2 + C n/2 2 , (4.1)
where we have set C = σ 2 β 2 e 2α -1 < 1. We deduce a simple algorithm to simulate the process σ between 0 and p. We fix the parameters N and n to get the desired approximation error and we compute the sequence (σ N,k ) 0 k n by applying iteratively the mapping T N . Note that it is sufficient to compute (σ N,n i ) i∈Z over the only dates i = -(n -k)N, . . . , p. In Figure . 1, we simulate both processes (σ i ) i∈Z and (r i ) i∈Z between the dates 0 and p = 300. (X n ) n is a stationary centered Gaussian process with covariation E[X i X i+j ] = λ 2 ln + (T /(j +1)) with λ 2 = 0.016 and T = 2000. The other parameters are chosen equal to α = 0.1, β = γ = 0.89, σ 2 = 0.025. We choose N = 300 and n = 10 in such a way that sup

k∈Z E[|σ k -σ N,n k | 2 ] ≈ 10 -12 .
5. Forecasting volatility 5.1. The general forecasting formula. We suppose that we observe (σ k , r k ) i 0 and we want to forecast σ i for i 1. For all i 1, we want to compute the best least squares linear predictor of σ i , i.e. the minimiser of:

inf α i (j),β i (k) E[(σ i -E[σ j ] - 0 j=-∞ α i (j)(σ j -E[σ j ]) - 0 k=-∞ β i (k)r k ) 2 ]
Figure 1. Example of simulation Thus, we introduce for all i, the coefficients (α i (j)) j 0 of the best least squares linear predictor of X i , namely the ones that minimise:

inf α i (j),β i (k) E[(X i - 0 j=-∞ α i (j)X j ].
We remind that one can find the (α i (j)) j 0 by solving the system:

C(|i -k|) = 0 j=-∞ α i (j)C(|j -k|), k 0 where C(j) = E[X i X i+j ] is the covariance function.
We can now state the main theorem of this section:

Theorem 5.1. Let us denote by σ i the best linear predictor of σ i . Then we have:

σ i = E[σ i ] + 0 j=-∞ α i (j)(σ j -E[σ j ]) + βσ 0 k=-∞ α i (k)r k (5.1)
where the α i (k) (independent of β) are given by α i (0) = -e -αi and for k -1:

α i (k) = ( 0 j=k+1 α i (j)e -α(j-k) ) -e -α(i-k) . (5.2)
In conclusion, we get prediction formulas that depend on C, α, β (a function and 2 parameters which have a clear signification). 5.2. An application: the long range correlation case. As an application of formula (5.1), we will compare the quality of the forecasts with and without Leverage for the multifractal model, i.e. E[X i X i+j ] = λ 2 ln + (T /(j + 1)). We study the SP500 index on the period 2002-2010. Our data set will consist of the daily returns r i and the following volatility proxy σ HL i on day i:

σ HL i = ln(H i /L i ),
where H i and L i are respectively the highest and the lowest value of the index on day i. In this context, we will use a filtering window of size L = 1000 (4 years) and use approximate formulas for the coefficients (α i (j)) -L+1 j 0 associated to the multifractal choice. More specifically, we choose the (α i (j)) -L+1 j 0 by discretization of continuous formulas first derived in [START_REF] Duchon | Forecasting volatility with the multifractal random walk model[END_REF] (i 1):

α i (j) = 1 π j j-1 √ i (L + i) √ -s (L + s)(t -s) ds, (5.3)
where L = 1000. To keep the paper self-contained, we give in the appendix a new derivation of formula (5.3). We will compare two set of values for α and β:

(1) The values in the absence of Leverage: β = 0.

(2) Fixing β = 5 (a typical empirical value: cf. With these two set of values, we will therefore compare the forecasts σ HL i derived from (5.1):

σ HL i =< σ HL > + 0 j=-999 α i (j)(σ HL j -< σ HL >) + β < (σ HL ) 2 > 0 k=-999 α i (k)r k ,
where the α i (j) are given by formula (5.3), the α i (j) by formula (5.2) and the < σ HL > and < (σ HL ) 2 > denote the empirical means:

< σ HL >= 1 1000 0 j=-999 σ HL j , < (σ HL ) 2 >= 1 1000 0 j=-999 (σ HL j ) 2 .
For the aforementionned set of values, we have plotted the renormalized empirical mean square error (MSE) as a function of the horizon (see figure 2 below where we set α = 1 Relaxtime ):

i → E[(σ HL i -σ HL i ) 2 ] E[σ HL i ]
.

We see that, for short horizons (depending on α), adding the Leverage effect improves the renormalized (MSE) whereas only certain values of α (around 0.01) improve the renormalized (MSE) for all horizons. This can be troublesome since the precise value of α that one would derive by fitting the Leverage correlation function seems to depend on the period (not to mention that the Leverage effect is very noisy and therefore the error bars on the estimation of α are huge!). 6. appendix 6.1. Proof of theorem 2.2. Let σ, γ, β be positive parameters and (X i ) i∈Z a square integrable stationary sequence. For 1 p < +∞, we introduce the metric space (E p , d p ) where:

E p = {(σ i ) i∈Z ; d p (0, (σ i )) < +∞, (σ i ) is F i -adapted and (σ i , X i , ǫ i ) i∈Z is stationary}
and the distance function is defined by:

d p ((σ i ), (σ ′ i )) = sup i∈Z E[(σ i -σ ′ i ) p ] 1/p
It is straightforward to check that (E p , d p ) is a complete metric space. We then introduce the mapping T defined on E p by:

(T σ) i = σ(γ + X i -β i-1 k=-∞ e -α(i-k) σ k ǫ k ). (6.1) For each element (σ i ) i∈Z , the sequence (σ i , X i , ǫ i ) i∈Z is stationary. It is plain to deduce that ((T σ) i , X i , ǫ i ) i∈Z is stationary. Provided that E(|X i | p ) + E(|ǫ i | p ) < +∞, we deduce that T maps E p into E p .
In the case p = 2, the mapping T is strictly contractive if and only if σ 2 β 2 e 2α -1 < 1. This can be seen by expanding:

E ((T σ) i -(T σ ′ ) i ) 2 = E σβ i-1 k=-∞ e -α(i-k) (σ k -σ ′ k )ǫ k 2 = σ 2 β 2 i-1 k=-∞ i-1 q=-∞ e -α(i-k) e -α(i-q) E (σ k -σ ′ k )(σ q -σ q )ǫ k ǫ q
The last double sum reduces to

σ 2 β 2 i-1 k=-∞ e -2α(i-k) E (σ k -σ ′ k ) 2 E[ǫ 2 k ]
because of the independence of (ǫ i ) i n from F n . It is then plain to deduce that

E ((T σ) i -(T σ ′ ) i ) 2 σ 2 β 2 e 2α -1 d 2 (σ, σ ′ ) 2 .
Therefore T is a contractive map if σ 2 β 2 e 2α -1 < 1. By the fixed point theorem, there is therefore a unique solution to the equation T (σ) = σ. Conversely, if there exists a square integrable non anticipative solution (σ, r) to (2.1) and (2.2), we can compute

E[σ 2 i ]
. By expanding the square as above, we obtain

E[σ 2 i ] =σ 2 (γ 2 + E[X 2 i ] + β 2 i-1 k=-∞ e -2α(i-k) E[σ 2 k ]) =σ 2 (γ 2 + E[X 2 i ] + β 2 e 2α -1 E[σ 2 i ]).
Therefore we get

E[σ 2 i ] 1 - σ 2 β 2 e 2α -1 = σ 2 (γ 2 + E[X 2 i ]) (6.2)
and we deduce σ 2 β 2 e 2α -1 < 1 (the process (X i ) i is assumed to be non trivial). Note that we have also proved that such a solution satisfies E[σ 2 i ] = σ 2 .

Remark: In the case when the processes (ǫ i ) i∈Z and (X i ) i∈Z are of p-th power integrable for p ∈ [1, +∞[, by using the Hölder inequality instead of expanding the square, we can show that for p ∈ [1, +∞[ the mapping T is strictly contractive for the metric d p if

σ p β p E[|ǫ 0 | p ] 1 e α -1 p < 1.
6.2. Proof of theorem 3.2. The proof is a simple adaptation of the proof in the discrete setting. For 1 p < +∞, we introduce the metric space (E, d p ) where E is the set of all random progressively measurable processes (σ t ) t∈R such that d p (0, (σ t )) < +∞ and (σ t , X t , (B t+h -B t )) t∈R is stationary for each h > 0 and where d p is the distance

d p ((σ t ), (σ ′ t )) = sup t∈R E[(σ t -σ ′ t ) p ] 1/p
(E, d p ) is a complete metric space. We define the mapping T c :

E p → E p by (T c σ) t = σ γ + X t -β t -∞
e -α(t-r) σ r dB r . (6.3)

In the case p = 2, it is straightforward to check that the mapping T c is strictly contractive if and only if σ 2 β 2 2α < 1. In that case, existence and uniqueness of a solution to (3.1) and (3.2) results from the Picard fixed point theorem. 6.3. Proof of theorem 4.1. We define the sequence (σ n ) n 0 of elements of E obtained by iterating the mapping T from 0, that is σ 0 = (0) i∈Z and ∀n 1, σ n = T (σ n-1 ) and we set C = σ 2 β 2 e 2α -1 . From classical arguments of fixed point theorem, we have the estimates

d 2 (σ, σ n ) 2 C n d 2 (σ, 0) 2 . (6.4)
Now we consider the mapping T N : (E, d 2 ) → (E, d 2 ) defined in section 4. Given a sequence u ∈ E, we can estimate the quantity d 2 (T u, T N u) as in the proof of Theorem 2.2 to obtain

d 2 (T u, T N u) 2 Ce -2αN d 2 (u, 0) 2 , (6.5) d 2 (T N u, T N u ′ ) 2 Cd 2 (u, u ′ ) 2 . (6.6)
Furthermore, by using (6.5) and (6.6), we have

d 2 (σ n , σ N,n ) =d 2 (T σ n-1 , T N σ N,n-1 ) d 2 (T σ n-1 , T N σ n-1 ) + d 2 (T N σ n-1 , T N σ N,n-1 ) C 1/2 e -αN d 2 (σ n-1 , 0) + C 1/2 d 2 (σ n-1 , σ N,n-1
), so that we get

d 2 (σ n , σ N,n ) n-1 k=0 C k/2 e -αN d 2 (σ n-1-k , 0). ( 6.7)
It remains to estimate d 2 (σ k , 0). We use the recursive relation

d 2 (σ n , 0) 2 =σ 2 (γ 2 + E[X 2 0 ]) + Cd 2 (σ n-1 , 0) 2 =σ 2 (1 -C) + Cd 2 (σ n-1 , 0) 2 ,
which is easily derived from the definition of T . We deduce

d 2 (σ k , 0) 2 σ 2 .
(6.8)

By gathering (6.4), (6.7) and (6.8), we complete the proof.

6.4. Proof of theorem 5.1. Since for all i, j E[X i r j ] = 0, we have:

σ i = E[σ i ] + σ 0 j=-∞ α i (j)X j -β 0 k=-∞ e -α(i-k) r k .
Since we have

X j = σ j σ -γ + β j-1 k=-∞ e -α(j-k) r k , we get: σ 0 j=-∞ α i (j)X j = 0 j=-∞ α i (j)(σ j -E[σ j ]) + βσ 0 j=-∞ α i (j)( j-1 k=-∞ e -α(j-k) r k ) = 0 j=-∞ α i (j)(σ j -E[σ j ]) + βσ -1 k=-∞ ( 0 j=k+1 α i (j)e -α(j-k) )r k
Finally, we get the prediction formula:

σ i = E[σ i ] + 0 j=-∞ α i (j)(σ j -E[σ j ]) + βσ 0 k=-∞ α i (k)r k
where α i (0) = -e -αi and for k -1, α i (k) = ( 0 j=k+1 α i (j)e -α(j-k) ) -e -α(i-k) . 6.5. The prediction formulas for fractional brownian motion and for the 1/f noise. We remind in this subsection the prediction formula for fractional gaussian noise dB H with Hurst index H ∈]0.5, 1[. Let S(R) denote the Schwartz space of C ∞ rapidly decreasing functions. We remind that the fractional gaussian noise is a centered Gaussian mesure in S ′ (R) (the space of tempered distributions) with covariance formally given by: E[dB H (s)dB H (t)] = H(2H -1) dsdt |t -s| 2(1-H) .

We restate the main theorem of [START_REF] Gripenberg | On the prediction of fractional Brownian motion[END_REF] (theorem 3.1): Theorem 6.1. We get the following prediction formulas (L¿0): where the kernel g H,L (s, t) is given by: g H,L (s, t) = sin(π(H -1 2 )) π t H-1 2 (L + t) H-1 2 (-s) H-1 2 (L + s) H-1 2 (t -s) .

To define the 1/f noise, one must introduce the space S 0 (R) of functions ϕ in S(R) such that R ϕ = 0. The 1/f -noise X is then the centered Gaussian measure in the quotient space S ′ (R)/R defined by (ϕ, ψ ∈ S 0 (R)):

E[ R ϕ(s)X s ds R ψ(t)X t dt] = R R ϕ(s)ψ(t) ln 1 |t -s| dsdt.
Since for all ϕ, ψ ∈ S 0 (R):

E[ R ϕ(s)dB H (s) R ψ(t)dB H (t)] 2(1 -H) → E[ R ϕ(s)X s ds R ψ(t)X t dt]
as H → 1, we can therefore recover the following prediction formula of [START_REF] Duchon | Forecasting volatility with the multifractal random walk model[END_REF] by letting H → 1: Theorem 6.2. We get the following prediction formula for the 1/f -noise X:

E[X t |(X s ) -L<s<0 ] = 0 -∞
g L (s, t)X s ds (6.9)

where the kernel g(s, t) is given by: g L (s, t) = 1 π √ t (L + t) √ -s (L + s)(t -s) .

We can finally recover formula (5.3) by discretization of formula (6.9).

[ 6 ] 1 Relaxtime

 61 ), we try different values of the Leverage correlation length: setting α = , we choose Relaxtime = 10, 30, 50, 200.

Figure 2 .

 2 Figure 2. Plot: empirical MSE for the SP500 index on the period 2002-2010.

E

  [dB H (t)|(dB H (s)) -L<s<0 ] = 0 -L g H,L (s, t)dB H (s)


On phase transition in classical fluid mixtures with surface adsorption


Simon Allier

Houari A Sahraoui

Salah Sadou
email: sadou@univ-ubs.fr

Stéphane Vaucher
email: vauchers@iro.umontreal.ca

Restructuring Object-Oriented Applications into Component-Oriented Applications by Using Consistency with Execution Traces

HAL is

Introduction

An intrinsic characteristic of software, addressing a real world activity, is the need to evolve in order to satisfy new requirements. Resulting from empirical studies, Lehman's first law states that software should evolve else it becomes, progressively, less satisfactory [START_REF] Lehman | Program evolution: Process of software change[END_REF]. Although old, this law has never been contradicted. The required reactivity (increasingly growing) of software applications, supports for business processes which are evolving more and more quickly, has even increased the scope of this law as the years go by. Maintenance is now, more than ever, an inescapable activity, the cost of which is ever increasing. Estimated at approximately 50 % to 60 % of the software total cost in the eighties and nineties [START_REF] Lientz | Problems in application software maintenance[END_REF][START_REF] Mckee | Maintenance as function of design[END_REF]; recent studies now evaluate this cost at being between 80 % and 90 % [START_REF] Erlikh | Leveraging legacy system dollars for e-business[END_REF][START_REF] Seacord | Modernizing legacy systems: Software technologies, engineering processes, and business practices[END_REF]. This high cost has undoubtedly been an effective catalyst for the emergence of new programming paradigms. Modular languages, then object-oriented languages, and more recently component-oriented programming, have always had as first justification, the significant increase in maintainability level. These new approaches can be used to build new applications. But what about legacy applications? In this case one can use the techniques of reverse engineering to transform the structure of the application, without changing its functionality, so that it conforms to the new paradigm.

In the past we presented a work that allows a company to organize its source code into reusable components [START_REF] Kadri | CBSE in Small and Medium-Sized Enterprise: Experience Report[END_REF]. In this work, identifying parts to put together in order to make a component was left in charge of engineers. What we propose is to automatically identify, in the case of object-oriented applications, classes to be grouped to form a component. Some approaches for identifying components [START_REF] Washizaki | A technique for automatic component extraction from object-oriented programs by refactoring[END_REF] or high-level architectures [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF][START_REF] Maqbool | Hierarchical clustering for software architecture recovery[END_REF] already exist. These works generally use the dependencies, between program elements, extracted by static analysis to identify highly cohesive and weakly coupled object-like structures. For components, the idea is to group together classes that contribute to the same functions. Generally, the restructuration of an application aims at improving its maintainability, including its evolvability. And often evolutions have a functional scope. But the classes used for building an application have in most cases a bigger scope than what is needed by the latter. So, unlike the approach based on static dependencies, we promote the dynamic dependencies for the aim of maintainability improvement. Thus, we clam that the most reliable way to determine which class contributes to which function is to execute the program with individual functions. Traces corresponding to execution scenarios could be analyzed to extract functional dependencies.

The rest of the paper is organized as follows: Section 2 describes the steps that constitute our approach. Scenarios execution and trace extraction are explained in Section 3. Section 4 gives details of our approach for identifying groups of classes that represent potential components. We provide three case studies in Section 5. Before concluding in Section 7, we describe some related works is in Section 6.

Approach Overview

We view a component as a group of classes collaborating to provide a system function. The interfaces provided and required by a component are the method calls respectively from and to classes belonging to other components.

So, we propose an approach for identifying components using traces obtained by executing scenarios corresponding to system use cases. The identification of candidate components consists of clustering the classes of the target system in such a way that classes in a group appear frequently together in the execution traces. In the same time, we also try to minimize the coupling between components.

Thus, the identification becomes an optimization problem where the goal is to identify groups of classes whose interactions are the most consistent with the execution traces. As this is a NP-hard graph partitioning problem, we use heuristic search to find a near-optimal solution.

Therefore, our approach is structured in three steps (see Figure 1):

1. Starting from a set of use case scenarios generate the execution traces. This step allows identifying dependencies between classes; 2. Produce a preliminary set of candidates using a global search. At this stage we use a genetic algorithm to produce this initial solution; 3. Refine the component candidates using a local search. We use a simulated annealing algorithm in order to achieve the local search. At the end of the identification process, groups that have a consistency with the traces below a predefined threshold are candidates to be packaged as components.

The execution traces are considered as the reference that guides the search for an identification solution. In the first step of the search, a population of potential groupings is created. Following an evolutionary algorithm, new groupings are derived that better match the interactions contained in the execution traces. This step serves primarily to find the region of the search space that has the highest potential to contain the best solution. Another step, in the form of a local search, allows to explore this region to find the near-optimal solution. This approach does not aim to identify reusable components. Indeed, component identification, in our case, is guided by the functional logic of the considered application. Restructuring the application using the identified components has the sole purpose of improving its maintainability. Thus, the goal of our identification process is not to fully re-architect a system into components. Our objective is rather to find groups of classes that could be packaged for reuse purpose.

It is true that in the case of a company, who work in a particular application domain, identified components can be considered reusable. But in this case, the execution traces must be obtained from several of its applications.

The extraction of the execution traces is described in Section 3. The two steps of the component identification are detailed in Section 4.

Execution Trace Generation

For the purpose of component identification, we are interested in finding dynamic relationships between classes. Hence, an execution trace (or call tree) is for us a directed tree T (V, E) where V is a set of nodes and E, a set of edges between nodes. Each node V i represents a class of the system (Cl i ). An edge V i , V j indicates that an object of class Cl i calls a method of an object of type Cl j . The root of T (V, E) corresponds to the entry point of the system or a thread run by the system. The execution traces are obtained by capturing the calls between instance of classes during an execution of a use case scenario. Every thread, created during the execution, produces an execution trace. In this initial version of the traces, nodes are labeled by both the actual types of the objects that are called and the methods called as showed in Figure 2(left). For example, when execution method m1, object B, instance of class B, called method m2 with as receiver object A, instance of class A. A, in turn, called m3, with as received D of class D.

In a second phase, nodes of the execution traces are relabeled by the classes corresponding to the called objects. The relabeling process is straightforward. Indeed as the call tree is generated dynamically, the concrete type of each object is recorded. This second phase produces traces where nodes are classes, and edges, dynamic relationships between classes. In the example illustrated in Figure 2 (right), the sequence <B.m1, A.m2, D.m3,> is replaced by <B, A, D>.

The identification of components from execution trace is relevant only if the execution traces cover all the functions of the system. Therefore, to extract the traces, we systematically apply all the recorded execution scenarios in the documentation.

Component Identification

As mentioned previously, the identification of components is modeled as a clustering problem. Indeed, the goal is to find the best partitioning P i (V ) = {C 1 , ..., C n } where V is the set of classes in the system and a C j is a a com-ponent candidate containing a subset of V . P i (V ) needs to satisfy the following completeness and consistency properties:

-0<j≤n C j = V (completeness) -∀j, k|j = k, C j ∩ C k = {∅} (consistency)
Exploring every possible clustering configuration cannot be done efficiently as we would need to consider every possible combination (NP-hard). We therefore propose using meta-heuristic search algorithms to find a near-optimal solutions. These algorithms are generic and are applied to a specific domain (in our case, graph partitioning) by defining the possible movements in the search space (transformations) and a fitness function that will be maximized.

In our approach, the search for a component-based architecture is implemented using a hybrid search [START_REF] Kelner | A hybrid optimization technique coupling an evolutionary and a local search algorithm[END_REF] which combines two different meta-heuristics: genetic algorithm (GA) and simulated annealing (SA). GA is a global search heuristic that applies changes to multiple solutions (called populations) and returns a solution that is globally near-optimal. This solution is then used as the initial solution of the SA algorithm, a local search algorithm. SA is called a local search algorithm because it explores the neighbourhood of a solution. Its output is our final solution. Both algorithms use the same solution representation and fitness function.

Solution and Space Representation

A solution for both algorithms can be any set of candidate components which respects the criteria of completeness and consistency. These candidates are represented by the set of classes they contain. In the example presented in Figure 2, the execution trace corresponds to a system with a set of classes A, B, C, D, E, F, H, I, J. One possible solution maybe: {{A, C, D}, {E, J}, {H, I, B, F }}. But there are many alternatives. In the following subsections we'll see how to choose the best alternatives.

Fitness Function

The fitness function used (Equation 1) evaluates the quality of a partition considering both the internal cohesion of components and the level of inter-component coupling of every component C (C can also be seen as a set consisting of classes that it gathers). The function takes as input A (for architecture), the set of component candidates proposed by the solution, and calculates the weighed average of the fitness of individual components (Cl being the set of classes in the system). The fitness of individual components depends mostly on their cohesion (Equation 2) unless the coupling level is too high in which case the fitness score is heavily penalized. The threshold cm used corresponds to the average coupling of all the classes in the system.

evalArch(A) = 1 |Cl| C∈A (evalComp(C) * |C|) (1)
evalComp(C) = evalCoh(C)/2 if evalCoupling(C) < cm evalCoh(C)/2 + 0.5 otherwise. (2)
Cohesion A good component should include classes that interact with one another to provide a specific set of functionalities; the strength of these interactions are what we call cohesion. The internal cohesion measure (Equation 3) evaluates how close are the different classes in the execution traces. It calculates the average distance between pairs of classes belonging to this component in all the traces. The distance between two classes a and b (Equation 4) is the average distance between instances of obj(a) of a and obj(b) of b.

To reduce the complexity of exploring the traces when calculating distances (number of edges), we use a constant d which indicates the maximum interesting distance between two classes. In other words, if we consider that the distance in the call graph, between two classes, that is acceptable to consider is for example d = 5, then we normalize the actual distance n using d (dividing it by 5). For distances less than 5, it will weigh less than 1 (acceptable) and for 5 and more is given a weight of 1 (worst penalty). The distance distM in(C, x, b) is the minimal distance between a node corresponding to an instance x of the class a and any instance of the class b in the execution trace where x appears. If no path with a length less than d is found, the distance is considered as d. all the node of the path are an instance of a class of C and such as the leng of the path is lower to

d, otherwise distM in(C, x, c) return d. evalCoh(C, d) =      1 |C| 2 -|C| x∈C y∈C,y =x dist(x, y, C) if |C| > 1 1 if |C| = 1 (3) dist(a, b, C) = 1 d|obj(a)| x∈obj(a) distM in(C, x, b) (4)
The number of classes into the same component increases with d. The choice of the best value of d can be done with small examples of applications from the domain and some experience in the domain.

Coupling One of the strengths of component-based development is that its components are loosely coupled and can be mixed and matched to build systems. The function evalCoupling(C) (Equation 5evaluates the level of coupling between components by counting the number of classes that are connected to a component (either calling or called). Classes that are part of the component are ignored.

evalCoupling(C) = | x∈C connected(x)| (5)
For the example presented in Figure 3, the solution S = {{A, B, C}, {D, F }, {E, H, I, J}} would produce the following metrics with d = 3: cm = 2.66 evalCoh({A, B, C}) = 0.47 evalCoh({D, F }) = 0.33 evalCoh({E, H, I, J}) = 0.7 evalCoupling({A, B, C}) = 2 evalCoupling({D, F }) = 2 evalCoupling({E, H, I, J}) = 4 evalArch(S) = ((0.47/2) * 3+(0.33/2) * 2+((0.7/2+0.5) * 4)) 9 = 0.49 In the following, we'll see how to determine the possible compositions.

Global Search

A genetic algorithm is a global meta-heuristics that emulates the concept of evolution. In order to find a solution from a population, it starts with a population (P 0 ) containing a set of solutions (called chromosomes) and simulates the passing of generations on this population. This initial population is a randomly generated. For every iteration of the algorithm, a new population (P i+1 ) is produced by selecting pairs of chromosomes ((c 1 , c 2 )) from P i and applying a crossover and/or a mutation transformation to these pairs with a certain probability. For this work, we systematically add the best chromosome of a generation to the next generation. The precise algorithm (Algorithm 1) uses three inputs: M axIter, M axN iter and M axSize. M axIter is the maximum number of generation for the evolution. M axN iter defines maximum number of generations where no improvement is accepted. Finally, M axSize defines the maximum size for a population.

Selection. The probability of selecting a chromosome c from the current population P depends on its quality with regards to the other chromosomes. This probability is given by the function:

P s(c, P ) = evalArch(c)
a∈P evalArch(a) This way of selecting components, called elitism, allows to give more chance to the fittest components to be selected.

Crossover. The "classic" crossover transformation consists of spliting two chromosomes c 1 and c 2 into two parts and merge the first part of c 1 (respectively the second part of c 1 ) with the second part of c 2 (respectively the first part of c 2 ). However, this crossover might generate a solution that does not respect the constraints of completeness and consistency. Indeed, some classes could exist in more than one component in the new generated chromosomes or do not exist at all. To preserve the two above-mentioned properties, we propose the following variation: Mutation. There are three type of mutation applicable to a chromosome:

-Split of a component in two components; -Merge of two components; -Move of a class from a component to another.

The type of the mutation is selected randomly as they are components involved in the mutation.

All three mutations produce solutions that preserve the properties of completeness and consistency. For example, c = {{A, C, D}, {E, J}, {B, F, H, I}} could be mutated in: c merge = {{A, C, D}, {B, E, F, H, I, J}} or into c move = {{A, C}, {D, E, J}, {B, F, H, I}}.

Local Search

GA can explore different solutions in a large search space to produces a solution that is globally near-optimal. This solution is then used by SA as a starting point for a fine grained exploration of its neighbourhood with the objective of refining it. The algorithm is presented Algorithm in 2. SA manipulate only one solution (s) at a time. At each iteration of the algorithm, this solution is compared to a neighbour (s neigh ) generated by a function N eigh(x). When s neigh is better than s as measured by the fitness function (evalArch), it replaces it. Otherwise, it can be accepted with a small probability which decreases as the algorithm progresses. This element of randomness is included to avoid falling into a local optimum.

Neighbour Function The neighbourhood function (N eigh(s)) uses the mutation of the genetic algorithm to produce a neighbour.

Case Study

In this section, we present and discuss the results obtained on three systems of different size (respectively 40, 73 and 221 classes).

System Descriptions

Our approach was evaluate on three systems. The first is an interpreter of the language Logo3 . It has a graphical interface which allows writing the code and a window which shows the result graphically. This programs contains 40 classes. Jeval4 is the second program. It is an expression interpreter. It contains 73 classes. Finally the last system is Lucene

Extraction of Traces

The extraction of execution traces was implemented using MuTT [START_REF] Liu | Feature location via information retrieval based filtering of a single scenario execution trace[END_REF]. MuTT is a Multi-Threaded Tracer built on the top of the Java Platform Debbuger Architecture. For a given program execution, MuTT generates an execution trace for each thread. For the three system, the extraction of the execution traces were generated as follows:

-Logo: The execution traces were obtained by executing different scenarios of use case. The definition of use cases and execution scenarios was easy, because, one of the authors of this paper was in the development team of the Logo interpreter. -Jeval, Lucene: The execution traces of Jeval and Lucene were derived by executing different scenarios of use cases. Scenarios was obtained by the test cases defined for these systems. We ensured that the test cases cover well the use-case scenarios.

Table 1 gives for every system the number of execution scenarios and the number of generated execution traces. There are more traces than scenarios because a trace is generated for each thread.

Result

The identification results for the three systems are presented in Table 2. For each system, it shows the number of identified components and the numbers of those who are related to the application and those who are not. A component is related if it contains the classes that provide a system function. It is considered as not related otherwise. Table 2. Results of component identification Logo Interpreter This system produces 5 components. Component Library implements the basic functionality of the language Logo (Math, String, ...). Display is composed of the classes responsible for the display of the instructions of the language Logo. Both components implement only one function. The third component Evaluator GUI provided two interrelated functions: the evaluation and the GUI for the result of the evaluation. The two other components do not contain functionally related classes. They both contain classes related to error management and other classes that plays the role of glue code between the three other components.

Jeval This system is partitioned into 5 components. Of these 5 components two represent respectively the library of the mathematical functions (sin, log, ...) and the library of the string functions. All the contained classes are related to the functions. The three other components contain classes necessary for respectively the parsing, the interpretation, and the mathematical operators evaluation. Depending on the viewpoint, these components could be merged.

Lucene From the 25 components, 16 are good and 9 bad. Four of the good components contain clear single functions. For example, QueryParser contain only classes responsible for the parsing of the search queries. The others 12 good components provide only one function, but with few classes missing. Here again, some could be merged if the goal is to obtained coarse-grained components. For example, the indexation function is split into subfunctions (5 components). Finally 9 identified components have no clear function.

As mentioned in section 2, the goal of our identification process is to find group of classes that could be packaged for reuse purpose. It is then normal that some of identified components are not considered as good. When putting a threshold on EvalComp, almost all the bad components will not be considered.

Components as a Behavioral Understanding Aid

Identified candidate components could be used to understand the behavior of a system. In the case of Logo Interpreter, when classes are grouped by their corresponding components in the execution traces, one can understand the behavior of the system. Indeed, the obtained nodes in the traces represent the system functions and, the links represent the function interrelations. Each execution scenario is then associated with a component interaction scenario. To illustrate the behavior understanding process, let us take the example of the following use case scenario:

Evaluation

Actor: Logo programmer Scenario: A1: Run the Logo interpreter A2: Write the following code in the editor window point 1 1 write "Hello" A3: Run the evaluation of the code from the editor window A4: Close the Logo interpreter

The execution trace corresponding to this scenario is shown in Figure 3. In this trace the component Evaluator GUI was manually divided into two components Evaluator and GUI. This trace contains two different phases : initial-ization and evaluation. In the first phase, when the Logo interpreter is run the GUI triggers the initialization of components Evaluator Library and Display (drawing standard output). In the second phase, the typed code is parsed and evaluated (Evaluator). During the evaluation, Evaluator calls twice Library: the first time for the function point and second time for the function write. The function point call the component Display to display the point at the coordinate (1,1), and the function call the component GUI to print in the standard output the text "Hello".

Discussion

The results of this case study are satisfactory. Indeed even if components identified are not all in the "related" category, the majority provide a unique feature and by spliting or merging the others, it is easy to obtain "related" component. Furthermore, this case study revealed a possible limitation of our approach. Indeed, our approach is designed to treat all the classes of the system as potential parts of components. It does not consider explicitly the case of glue-code classes. Detecting such classes and excluding them from the partitioning will certainly improve the identification results.

Related Work

The work proposed in this paper crosscuts three research areas: architecture recovery/remodularization, legacy software re-engineering, and feature location.

Different approaches have been proposed to recover architectures from an object-oriented program. The Bunch algorithm [START_REF] Mitchell | On the evaluation of the bunch search-based software modularization algorithm[END_REF] extracts the high-level architecture by clustering modules (files in C or class in C++ or Java) into subsystems based on module dependencies. The clustering is done using heuristicsearch algorithms. In [START_REF] Medvidovic | Using software evolution to focus architectural recovery[END_REF], Medvidovic and Jakobac proposed the Focus approach whose goal is to extract logical architectures by observing the evolution of the considered systems. The approach identifies what the authors call processing and data components from reverse engineered class diagrams. Closer to our work, the ROMANTIC approach [START_REF] Chardigny | Quality-driven extraction of a component-based aachitecture from an object-oriented system[END_REF] extracts component-based architectures using a variant of the simulated annealing algorithm. In addition to dependencies between classes, other information sources are considered to help identifying functional components rather than logical sub-systems. Such sources include documentation and quality requirements. In [START_REF] Washizaki | A technique for automatic component extraction from object-oriented programs by refactoring[END_REF], Washizaki and Fukazawa concentrate on the extraction of components by refactoring Java programs. This extraction is based on the class relation graphs. In the above-mentioned work, the component extraction process uses dependencies between classes/modules that are extracted using static analysis. Dependencies are not related to particular functions of the considered system which makes it difficult to relate identified components to specific functions. Thus, the components identified have a general scope and are not dedicated to the application. Whereas in the case of our approach, identification is guided by the context of the application, which will facilitate its maintenance.

Our approach use heuristic-search methods [START_REF] Kelner | A hybrid optimization technique coupling an evolutionary and a local search algorithm[END_REF], genetic algorithm and simulated annealing. Search-based methods is widely applied to solve problems similar to ours. For example, Seng [START_REF] Seng | Search-based determination of refactorings for improving the class structure of object-oriented systems[END_REF] improves the design of an 00 code with a fitness function that quantify the improvements in design metrics. To this end, a genetic algorithm is used. In [START_REF] Kessentini | Model transformation as an optimization problem[END_REF], Kessentini uses meta-heuristiques to transform models by examples. More close to our work, [START_REF] Chardigny | Quality-driven extraction of a component-based aachitecture from an object-oriented system[END_REF] use a variant of the simulated annealing for extract component-based architectures.

Feature location is probably the problem that is closest to the one addressed in this paper. Many research contributions proposed solutions that are based on dynamic analysis [START_REF] Greevy | Correlating features and code using a compact twosided trace analysis approach[END_REF] or combinations of static and dynamic analyses [START_REF] Eisenbarth | Locating features in source code[END_REF][START_REF] Liu | Feature location via information retrieval based filtering of a single scenario execution trace[END_REF][START_REF] Poshyvanyk | Feature location using probabilistic ranking of methods based on execution scenarios and information retrieval[END_REF]. In general, static analysis uses call graphs and/or keyword querying. From the other hand, dynamic analysis consists in executing scenarios associated with features, collecting the correspondent computational units (PU), generally methods or classes, and deciding which PU belongs to which feature. The decision can be made by metrics [START_REF] Greevy | Correlating features and code using a compact twosided trace analysis approach[END_REF], FCA [START_REF] Eisenbarth | Locating features in source code[END_REF], or a probabilistic scoring [START_REF] Poshyvanyk | Feature location using probabilistic ranking of methods based on execution scenarios and information retrieval[END_REF]. Finally, sometimes, static analysis is used to enrich the dynamic analysis. Both analyses can also be performed independently and their results combined using a voting/weighting function. The combination of static and dynamic analyses is also used in a semiautomatic process where visualizations are proposed to experts to make their decisions [START_REF] Bohnet | Visual exploration of function call graphs for feature location in complex software systems[END_REF]. Like in our case, the feature location approaches use dynamic analysis and try to associate program units to scenarios. The difference, however, is that the problem of locating features and identifying components are different in objectives and nature. In the first case, the goal is to determine code blocks, methods, or classes that are involved in a particular feature represented by a set of scenarios. For component identification, a scenario may involve many features (data acquisition, processing, data store, and results displaying). The association between feature and scenario is not a one-to-one relation. Moreover, the execution of a feature may necessitate the execution of many other features, which makes it difficult to draw the boundaries. For this reason, we view the component execution as sequences of interactions between classes in an integrated dynamic call graph.

Conclusion

Our main objective in restructuring an object-oriented application into a component-oriented application is the improvement of its evolvability. Thus, unlike other existing approaches where dependencies between classes are derived by static analysis, we used, for our component identification approach, method call trees obtained by executing use case scenarios on the application. This guarantees that only functional dependencies are considered in the component identification. Indeed, application's evolutions have most often a functional scope. Moreover, the execution traces, obtained thanks to the use cases, limit the analysis of dependency only to the space covered by the application. While the classes, which are often generic, cover a wider space.

In the past we presented a work that allows a company to organize its source code into reusable components [START_REF] Kadri | CBSE in Small and Medium-Sized Enterprise: Experience Report[END_REF]. This consisted in a reorganization of the development life-cycle and in the use of the UML2 component model in order to wrap a code corresponding to a component. Although in the case of our experimentation with our industrial partner, the engineers know very well the existing code, identification remains a tedious job. Moreover, in several cases of incorrect component identification, the cause was a reflex related to their long experience with the object-oriented approach.

The work presented here is complementary to that presented above. It is an aid in identifying components by grouping classes. As our approach relies on execution traces of an application, the proposed grouping is necessarily adapted to this application. What goes in the direction of improving the maintainability of the application. In the case of our industrial partner, applications are built to the same scope (Geographical information systems) using their library of classes. To build their library of components, we should grouping classes by using traces from all their applications. After that, we plan to use another work, that we have already made, in order to automatically select components [START_REF] George | A methodological approach for selecting components in development and evolution process[END_REF].

Our approach comes just before building the components. It only proposes the classes that must go together. Thus, its use is possible with any model of components.

Fig. 1 .

 1 Fig. 1. A three-step process for component identification

2 2 Piter+1 1 = 2 =

 2212 Algorithm: genetic(M axIter, M axN iter, M axSize) let iter := 0; niter := 0 create a initial population P0 let Best := minc∈P evalArch(c) while (iter < M axIter) and (niter < M axN iter) do eval Piter let Piter+1 := {∅} while Piter+1 < M axSize do Select c1, c2 ∈ Piter Crossover c1, c2 with probabilitypc to c ′ 1 , c ′ Mutate c ′ 1 , c ′ 2 with probability pm to c ′′ 1 , c ′′ := Piter+1 ∪ {c ′′ 1 , c ′′ 2 } end let BestLocal = minc∈P iter evalArch(c) Piter+1 := Piter+1 ∪ {BestLocal} if evalArch(BestLocal) < evalArch(Best) then Best := BestLocal niter := 0 end iter + +; niter + + end return Best Algorithm 1: Genetic algorithm -Divide the chromosome c 1 (respectively c 2 ) into two parts c 11 and c 12 (respectively c 21 and c 22 ), each containing a subset of components. -Create a chromosome c ′ 1 by insert c 11 between c 21 and c 22 (respectively c ′ 2 by inserting c 21 between c 11 and c 12 )). -Delete in c 21 and c 22 all the classes appearing in c 11 (respectively in c 11 and c 12 all the classes appearing in c 22 ). For example, the chromosomes: c 1 = {{A, C, I}, {E, J}, {D, H, B, F }} c 2 = {{A, H}, {B, C, D, E}, {F }, {I, J}} partitioned into: {{A,C,I}} and {{E, J}, {D, H, B, F }} for c 1 , {{A,H}, {B,C,D,E} and {F }, {I, J}} for c 2 . produces the two chromosomes: c ′ {{H}, {B, D, E}, {A,C,I}, {F }, {J}}, and c ′ {{I}, {A,H}, {B,C,D,E}{J}, {F }}.

Fig. 3 .

 3 Fig. 3. Execution trace

Table 1 .

 1 5 , a high-performance, full-featured text search engine. Lucene contains 221 classes. Capture of the executions traces

		number of executions trace
	Logo	8	19
	Jeval	9	9
	Lucene	19	59

http://naitan.free.fr/logo/

http://jeval.sourceforge.net

http://lucene.apache.org/


Sablayrolles B Azevedo-Coste

C Demar

Torque Prediction Based on Evoked EMG in Fatiguing Muscle Toward Advanced Drop Foot Correction

Keywords: drop foot correction, muscle fatigue, stimulus evoked EMG, electrical stimulation

Electrical stimulation (ES) has been applied since 1961 for the correction of hemiplegic drop foot. One main drawback of the technique is the occurrence of early fatigue. Several studies attempted to propose solutions to decrease the fatigue, however no consensus was reached to date. This work aims to predict ankle torque using stimulus evoked EMG (eEMG) during different muscle fatigue states. Five healthy subjects participated in our study. Conventional stimulation for drop foot correction was applied by surface stimulation in sitting position. The results showed that no tendency between the changes in eEMG and torque with muscle fatigue was found among all subjects. However, the relationship variation in one subject gradually changed depending on muscle states. In this work, we carried out the torque prediction with an adapted parameters model according to muscle fatigue state by reidentification using the latest measurement. The prediction was improved with 21%~90.9% comparing to the fixed parameters model. The results revealed a promising approach to use evoked EMG for fatigue compensation in the application of drop foot correction.

Introduction

Hemiplegic drop foot (DF) is a condition resulting from stroke where the subject cannot lift his or her foot or dorsiflex [START_REF] O'keeffe | The Development of a Potential Optimized Stimulation Intensity Envelope for Drop Foot Applications[END_REF]. The lack of dorsiflexion results in a dragging of the leg and as a result has a significant impact on the person's gait [START_REF] O'keeffe | The Development of a Potential Optimized Stimulation Intensity Envelope for Drop Foot Applications[END_REF]. Electrical stimulation (ES) is one of the existing solutions which are used to correct drop foot. Many researchers paid most attention to get optimal stimulation intensity profiles for the affected leg. An optimised stimulation envelope to reproduce the EMG pattern observed in the tibialis anterior (TA) during healthy gait was proposed in [START_REF] O'keeffe | The Development of a Potential Optimized Stimulation Intensity Envelope for Drop Foot Applications[END_REF]. The changes in TA activity with walking speed was further considered to enable adaptive FES intensity envelope in [START_REF] Byrne | Effect of Walking Speed Changes on Tibialis Anterior EMG during Healthy Gait for FES Envelope Design in Drop Foot Correction[END_REF]. A method providing with precise timing adaptation of the ES pattern to ensure a good coordination of the healthy and affected legs was proposed in [START_REF] Azevedo-Coste | Towards Continuous Online Adaption of Electrical Stimulation Sequences for Drop Foot Correction[END_REF]. In practice, the clinical implementation of ES for drop foot correction in hemiplegia is a challenging task. One of the challenges is the rapid onset of fatigue on the stimulated muscle. Thus, it is required to adjust stimulation parameters to provide optimal and tolerable treatment with minimal muscle fatigue. In previous researches, optimal stimulation patterns or modulation of stimulation parameters were proposed to reduce fatigue. However, no consensus was reached to date. Moreover, fatigue phenomenon is unavoidable with repetitive or prolonged stimulation and it usually occurs fast with ES on paralyzed muscles. Therefore, it is essential to predict force generation for precise ES closed loop control when the stimulated muscle becomes fatigued.

The aim of this work is to predict ankle torque preliminarily in healthy subjects for the application of ES induced ankle dorsiflexion. In our study, we propose to identify the eEMG-torque model adaptively to improve torque prediction.

Material and Methods

Experimental Setup

Five healthy subjects participated in this study. A commercial stimulator (Prostim) was used to induce sufficient dorsiflexion. The active (cathode) stimulating electrode was placed over the common peroneal nerve and the indifferent (anode) was placed over the TA muscle. The isometric ankle torque was measured by a calibrated dynamometer (Biodex), interfaced with an acquisition system (Biopac MP100). The subjects were seated on the chair with their right ankle at 90 o , while the foot was strapped on the pedal. Evoked EMG activity of TA was collected using bipolar surface electrodes. The electrodes were positioned on the TA muscle along muscle fiber direction with 20mm interelectrode spacing. The reference electrode was placed on the patella. The EMG signal was amplified (gain 1,000) and sampled at 4KHz.

Experimental Protocol

For each subject, the first step of the experiment is to find the optimal stimulation location, that is, inducing pure dorsiflexion but tolerable and painless. The second step is to find supramaximal stimulation amplitude. In order to imitate hemiplegic walking situation, sequences of 2s stimulation and 2s rest were applied during 30mn. Stimulation consists of a trapezoidal envelope with pulse width modulation. Each 2s stimulation train consists of 0.4s ramp-up, 1.2s plateau and 0.4s ramp-down. The maximum pulse width and stimulation frequency were fixed at 350 s μ and 40Hz for all subjects.

Signal Processing

The stimulation artifacts were removed by means of blanking window to extract muscle response (Mwave) signal. The measured torque was offset with respect to the baseline of the torque measurement without stimulation. The torque data was lowpass filtered (6 th order, cutoff frequency 100Hz). The measured eEMG data was lowpass filtered (6 th order, cutoff frequency 300Hz). The eEMG data was divided into epochs with each epoch containing one Mwave. The mean absolute value (MAV) of eEMG was calculated every 9 Mwave through moving average. The mean torque during the same time window was simultaneously computed.

The measurements of torque and eEMG were used to identify the parameters of mathematical models for ankle torque prediction. A Hammerstein model was adopted to represent the contraction dynamics as proposed in [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically stimulated muscle[END_REF]. The eEMG-torque model is shown as follows:

Where y(t) is the torque output at time t, x(t) is the input, in this study, it is the MAV value of eEMG. The term c o is used to fit any offset of the output torque. We choose n=3 as in [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically stimulated muscle[END_REF]. Different model order (l,m) was assumed to identify and validate the model. Finally we choose model order l=3, m=4 to minimize the prediction error. A recursive least squares method was used to fit these parameters.

For model identification, past measured torque is used as the past torque y(t-i). For torque prediction, past predicted torque is used as the past torque. In this case, model output

) (t y p )
based on the identified parameters can be computed as follows: When force measurement is not available, this approach makes it possible to use evoked EMG as a synthetic force sensor. In this case, we initialized the predicted torque at zero when no stimulation was delivered to the muscle. Mean-squared prediction error (MSE) is used as prediction performance index.

Results

For each subject, in order to observe the global change tendency in torque and eEMG, the peak torque and peak MAV during every 2s stimulation train were calculated and then normalized against the maximum values during the whole stimulation duration. We did not find consistent relationship between the variation in eEMG and torque associated with muscle fatigue among these subjects. Fig. 1 showed the changes in torque and MAV of eEMG with fatigue in subject1. The initial 2mn stimulation resulted in significant potentiation of muscle activity, whereas the torque gradually decreased. From 1.4mn to 4.1mn, the torque maintained at the same level. However, the MAV underwent three different states: increase, transitional remain and finally decrease. We also found time-varying properties in the variation of eEMG and torque in all other subjects. This kind of dissociation due to muscle fatigue condition was also found in [START_REF] Thomas | EMG Changes in Human Thenar Motor Units With Force Potentiation and Fatigue[END_REF]. Therefore, it is difficult to obtain precise torque prediction using a fixed eEMG-torque model especially in prolonged stimulation. In order to validate the model, we adopted two methods: 1) fixed parameters model as proposed in [START_REF] Erfanian | Using Evoked EMG as a Synthetic Force Sensor of Isometric Electrically stimulated muscle[END_REF] and 2) our method based on adapted parameters model. For each subject, we divided the eEMG data into sets which contain 6 trains (33s). Then we subdivided every set into two epochs with each epoch containing 3 trains (16.5s). In the first approach, we applied the initial 16.5s data in first set to identify model and fixed the model parameters to predict all the following ankle torque. In the second approach, we applied the data during the first epoch (16.5s) in each set to identify model. Then the obtained parameters were used to predict the torque just during the adjacent epoch (16.5s) in the same set. The process was repeated for every set during the whole stimulation duration. Summary of the prediction performance is shown in TABLE I. Test1~Test4 were selected at different time instances to have different muscle states in terms of muscle fatigue. We can find that the prediction is less precise or difficult with the fixed model (e.g. test2, test3 in subject1) during long-term repetitive stimulation. However, when the model was reidentified with the latest data as our approach did, the prediction errors were greatly reduced as shown in TABLE I.

In figure 2, we selected three sets of data (Test1, Test2 and Test 4) at different muscle fatigue states in subject2 to compare the prediction performance of both approaches. MSE1 and MSE2 respectively indicate the prediction error with fixed and adapted model. 30mn stimulation resulted in torque loss from (a) 100%, to (b) 60%, until (c) 20% mainly due to muscle fatigue. The prediction with fixed model became less precise over time. Nevertheless, the torque prediction was greatly improved with adapted model.

Discussion and Conclusions

Our study focuses on the prediction of ankle torque based on evoked EMG in the context of ES drop foot correction during different muscle fatigue states. 5 healthy subjects participated in the experiments. During 30mn stimulation applied in isometric conditions, the torque declined in all subjects. The change of eEMG did not represent tendency among them. However, eEMG-torque relationship was gradually varied according to muscle fatigue condition in the same subject. Despite this kind of dissociation between eEMG and torque, we were able to apply our adapted eEMG-torque model to highly improve prediction performance comparing to fixed model which was commonly used in previous researches. Thereby, the effect of adapted identification should be significant and this approach has great potential for contribution of ES closed loop control with fatigue consideration. The future work includes to validate the results in hemiplegic subjects, to extend the present method to on-line automatic identification and adaptive control of ES for fatigue compensation towards advanced drop foot correction.

Fig. 1 :

 1 Fig. 1: The transition of the torque response and MAV of eEMG by stimulation in subject1. The results are normalized by the maximum values for the whole test.

Fig. 2 :

 2 Fig. 2: The measured torque and predicted torque based on evoked EMG in different muscle fatigue state in subject2. The measured torque (blue), the predicted torque with fixed model (green dotted) and the predicted torque with adapted model (red dashdotted) are shown.