F Cottet
email: cottet@ensma.fr

L David
email: david@ensma.fr

IMPROVING TEMPORAL BEHAVIOR WITH GRAPHICAL METHOD IN REAL-TIME SYSTEMS

Keywords: validation of real-time systems, timing analysis, fault-tolerant software, scheduling algorithms

The purpose of this work is to provide to the designer a graphical method for taking design decisions for enhancing the temporal performance of the application in terms of fault tolerance, task response time, etc. We have introduced a graphical representation of the schedulability analysis based on a classical temporal model in a given scheduling environment. This methodology is extremely useful for supporting both control applications and multimedia systems, in which the execution times or rates of some computational activities can vary significantly, and for which the timing fault tolerance must be evaluated in order to prevent failures.

INTRODUCTION

A real-time system is, by definition, a system that must satisfy explicit response-time constraints, or risk severe consequences including failure [START_REF] Laplante | Real-Time Systems Design and Analysis -an Engineer's Hand Book[END_REF]. Since most of the mission-critical and human-critical systems operate in real-time, these also require a high degree of dependability. It can then be sustained that two of the most important characteristics of real-time systems are the capability to fulfill strict timing requirements and to guard against imperfect execution environments.

In real-time systems, the time is a main point, unlike the definition for non-real-time systems, where time is only a performance matter. So an appropriate approach is a temporal analysis and validation of the system before its execution in order to guarantee the timing constraints of the application. One method, increasingly used, is based on schedulability analysis. Various methodologies are founded either on specific real time languages [START_REF] Stoyenko | The evolution and State-of-the-Art of Real-Time Languages[END_REF], or on a temporal modeling of software written with a classical high level language associated with a real-time kernel (Cottet et al., 1994;[START_REF] Zhou | rate-Monotonic Scheduling in Presence of Timing Unpredictability[END_REF], or on formal methods that lead to exhaustive computation of the scheduled task execution sequences [START_REF] Geniet | Exhaustive Computation of the Scheduled Task Execution Sequences of a Real-Time Application[END_REF]. This operational validation by schedulability analysis forces the mathematical characterization of the task set according to a precise temporal model. In the simplest widely accepted model, a task Ti, assumed to be periodic and independent, is represented by the quadruple (ri, Ci, Di, Pi) where ri is the first request of the execution of the task, Ci is the computation time, Di is the deadline and Pi is the period. We have the basic obvious conditions (1) (Ci ≤ Di ≤ Pi ) where Di ≤ Pi means that the task deadline is lower than the next execution request (the task must be completed before a new invocation) and Ci ≤ Di the necessary condition for the execution.

Three of these parameters (ri, Di, Pi) are elicited from the timing requirements giving the temporal characteristics of the input/output operations and the process behavior. The temporal parameter Ci corresponds to the execution time of each task in the intended hardware and software execution environment. It can be evaluated by different ways : either from an analysis of the program code [START_REF] Puschner | A Tool for the Computation of Worst Case Tool Execution Times[END_REF] or from a direct measurement of execution time (Babau et al., 1995). It is usually fixed and depends only on the algorithmic implementation in a given hardware context. It is possible to estimate a variation range from a minimum run-time Ci,0 to the worst case execution IFAC Workshop On Real Time Programming, Palma de Mallorca, Spain, April, 2000, pp. 79-84 time Ci,m of a task. The others parameters can be adjusted by the designer to obtain a system that matches as closely as possible the timing constraints of the application. Particularly, this adjustment of the deadline and the period of the task is extremely important because, in relation to the formal scheduling algorithms such as Rate Monotonic and Earliest Deadline [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF], it corresponds to a priority mapping that takes into account the importance of tasks.

The purpose of this work is to provide to the designer a graphical method of task temporal parameter adjustment in order to greatly improve the temporal behavior of a hard real-time system. This approach is particularly useful to enhance the performance of the application in terms of quality of temporal service (task response time, processor idle time, number of processor context switches, etc.). Moreover this graphical method can be used to provide tolerance of timing faults caused by the effect of the parameter variation as the period jitters or the computational time uncertainty. It provides an approach to reliable real-time systems. The main investigations for improving fault tolerance of hard real-time systems were realized by modifying task model [START_REF] Kim | Reliability Modeling of Hard Real-Time Systems[END_REF]. The paper is organized as follows. In section 2, we present the graphical model of a task. In particular, we choose to describe more precisely the plane "deadline versus period". In section 3, we introduce the schedulability condition in the graphical representation. A simple example ends this section. Section 4 presents the graphical analysis of the plane "execution duration versus period". A speediness increase of this methodology and a real industrial example are presented in the section 5. Conclusion are given in section 6.

GRAPHICAL MODELING OF TASKS

Scope and assumptions

This work is focused on one-processor and preemptive real-time systems. Concerning the tasks, we simply identify two different kinds of activation signals. The first concerns the periodic signals, characterized by a period P, which activate hardware periodic tasks. The second class contains the nonperiodic signals that occur at arbitrary times without a fixed time interval. In this situation, we need a classical specification that is given by the minimum interval dmin between two occurrences of such signals For this latter case, we have adopted the simplest approach very practical in real-world applications : the immediate periodic server (the server task services any pending non-periodic requests, if no non-periodic requests are pending, the server suspends itself until its next period and the released time is instead used for periodic tasks) (Homayoun and Ramanritham, 1994;[START_REF] Spuri | Scheduling aperiodic tasks in dynamic priority systems[END_REF][START_REF] Buttazzo | Hard Real-Time Computing Systems, Predictable Scheduling Algorithms and Applications[END_REF]. So from now on, we will consider that a periodic server is created and dedicated to each non-periodic request. Therefore we assume that a real-time software is depicted by periodic task set.

A three-dimensional general model of a task

Assuming that each task has a known initial invocation ri, the key point is to skillfully determine the other temporal parameters (Ci, Di, Pi) according to the non-functional requirements such as the timing properties of the application. First, the parameters of the task temporal model (Ci, Di, Pi), as previously presented, must satisfy the relations (1). Therefore, in a set of axes (run-time, deadline, period), the three-dimensional graphical representation of the volume corresponding to the allowed values of task parameters has a trihedral shape. If we draw the projections of this sphere on the different planes, we obtain the areas of the allowed values of task attributes : plane projection (period versus run-time), front elevation (deadline versus period) and side elevation (deadline versus run-time). All theses surfaces come in the form of an opened wedge. The vertices of these acute angle wedges correspond to the equality of the parameters: Ci = Ri = Pi = Ci,0. So, if we assume that the computational time varies in the range [Ci,0, Ci,m], then, in the deadline versus period plane for example, the allowed zone changes as described on the figure 1. The more the run-time increases, the more the possible surface of task parameter values shrinks.

We confine our attention here to two of these projections that seem quite interesting for our purpose of temporal analysis of task behavior : deadline-period and period-execution time coordinate planes. We first assume that each task has a known worst case upper bound to the execution time that can be computed (Ci=Ci,m). So the free other parameters (Di, Pi) can be studied according to the timing requirements of the task itself, derived from the global timing properties of the application. The effect of relaxing the first assumption (Ci=Ci,m=constant) will be discussed in section 4.

Graphical temporal modeling of periodic tasks

Up to now, the task was considered as a theoretical model of periodic task without taking into account some timing requirements of the application. From these temporal constraints related to the controlled process, different temporal characteristics of a task can be deduced and imposed ; for example, we can have the following requirements :

D : Deadline

P : Period

Run-ability constraint (line D=P)

C i,0 C i,0 C i,m C i
• a maximum period Pi,max (Pi≤Pi,max), corresponding, for instance, to a sampling rate limit, • a minimum period Pi,min (Pi≥Pi,min), corresponding to a useless repetition rate (for instance the period of a polling task with a value much greater than the physical signal change rate), • a maximum deadline Di,max (Di≤Di,max), in order to get a better response time.

The graphical representation of figure 2 shows the previously determined area of the allowed values of task attributes located inside the shaded wedge, named zone A. On the other hand we obtain a more restricted surface (zone B), included in the first one, that corresponds to the set of potential couples (Di, Pi) of a given periodic task satisfying its timing constraints.

-Temporal modeling of server tasks

A non-periodic event is assumed to be specified by at least a minimum time interval between two invocations, d min , and if necessary a maximum deadline, Dev. So we have to find the temporal parameters (Dsi, Psi) of the periodic server that offers a correct service to any of those non periodic invocation signals, that is to say : to respond with a delay less than its deadline Dev. As before, we assume that the computation time Csi allocated to the immediate server is known. Given that there is of course no synchronization between the periodic server activation and the occurrences of the non periodic events, we are going to consider the worst case. Suppose that the non periodic event occurs just after the release time of the periodic server and that it suspends immediately since there is no pending nonperiodic signal. Therefore the non periodic request event will be served during the next period of the server but before the server deadline attached to this new period. Due to that, the temporal parameters must meet the following condition :

Psi + Dsi ≤ Dev ≤ dmin (2)
This condition, involving Psi and Dsi, offers a certain choice of the periodic server parameters : either a long period server with a short deadline that strongly constraints the schedulability of the task set or a short period with a deadline equal to the period that is prejudicial to the processor load.

-SCHEDULABILITY HANDLING IN THIS GRAPHICAL MODELING

The temporal model representations so far show the allowed main temporal parameters of the task with the effect of the timing requirements, but without taking into account the influence of other tasks, or, in others words, the concurrent execution context. The aim of this work is to provide a methodology of task set analysis for determining the well adapted temporal parameters. The best fitting parameters are considered with respect to the global timing requirements of the application. This property requires a framework that incorporates the execution environment, especially the scheduling algorithm.

We have first chosen a static scheduling algorithm (the priorities are not changed once they are assigned to the tasks) namely the Rate Monotonic (RM) algorithm. This scheduling algorithm is a simple rule that assigns priorities to tasks according to their request rates. Specifically, tasks with shorter periods will have higher priorities.

Let us consider a real-time system with k tasks (basic periodic task and/or server task). The worstcase computation time Ci of each task Ti is assumed to be known in a deterministic manner. Our approach to fixing temporal task parameters of the task set works according to an iterative process. At a given step, we assume that the parameters of k-1 tasks have been determined, that is to say, their point (Di, Pi) or (Dsi, Psi) has been chosen in the allowed areas. From these assumptions, we can build the graphical representation (Dk, Pk) of the task k with its proper restricted surface and the boundaries representing the interference of the concurrent execution context caused by the k-1 existing tasks. This methodology can be seen according to two different points of view. Either, in a design process, the k-th task must be added to the others or, in a redesign or analysis process, the k-th task is one among the others whose temporal characteristics one wants to improve.

The first obvious bound is given by the maximum processor utilization factor, that is to say 1 in a oneprocessor context. Therefore, Pk has a minimum value Plim which is given by :

P k P lim C k (1 C i P i i 1 k 1 (3)
But the main limitation that we have to take into account is the schedulability constraint. For a RM scheduling and a task set where each task must be completed before the new invocation, a whole task set is schedulable if [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real-time environment[END_REF] :

C i P i i=1 k-1 C k P k k 2 1/k 1 (4)
In this particular assumption (run-ability constraints for all the tasks) the preceding condition gives us a minimum period PLL,lim located on the "D = P" line.

It should be noted that this condition is sufficient but not necessary. It should be pointed out that these different limitations are often very pessimistic. To improve efficiency of this computer aided design methodology, we have to do a simulation of the RM scheduling of the task set. Then the designer can choose the appropriate pair of temporal parameters (see figure 3). This simulation must cover scheduling for a length of time equal to the major cycle corresponding to the least common multiple of the task periods LCM{Pi}.

In the case of a server task (example of the figure 3), a very important limit can be defined : the maximum response time TR of the task that corresponds to the sum of the deadline and the period. In that case also, this limit is an upper bound, the server task can be replied more quickly to the aperiodic interruption in a specific execution sequence. So, for a given scheduling algorithm, the direct computation of valid sequences permits us to draw the real and complete limits of the allowed zone for a given task in a concurrent scheduling environment. This final zone (C), that is included in the zone B (previously inside the zone A), can be much more restricted than the previous surfaces. The borders of this surface are constituted by of different lines : the D=P line where we find the analytical results for the period, the specification upper value of the deadline, the maximum response line, the minimum deadline that can be irregularly shaped. The final zone C emphasize the capability of the system to tolerate timing fault related to the deadline and the period of this particular task.

First representing the results of an simple example in the D-P plane is slight different from the preceding theoretical graphs because the D and P variables are integers and therefore the D-P plane is discrete. Let the test example be a periodic task set with two tasks characterized by the following temporal parameters : (0, 2, 8, 10) and (0, 2, 8, 8). Assume that we want to add a periodic server task created in order to deal with a non-periodic event, characterized by a needed computation time C3 of 3 and a minimum occurrence interval dmin of 15. We can calculate the minimum period corresponding to the fully used processor according to the equation ( 3) : Plim = 6. To obtain the complete graph, the RM scheduling has been applied for the predetermined range of P3. The figure 4 shows the final results with the allowed zone of temporal parameters of this task.

The designer can chose among all the points the possible values of D and P. There is no better choice among all these possible parameters. Different criteria can be used to select the task temporal characteristics such as minimizing the processor load : point A(10, 5). Another way is to optimize two criteria like processor utilization and response time. In this case, the point, referenced B(8, 7) on the figure 4, seems to be a very interesting compromise between the highest available period and a correct response time.

In all the above execution simulations, we consider a RM scheduling. Before modifying the temporal parameters of a task, a designer could also want to analyze how a more efficient scheduling policy can widen the valid zone. So we now turn our attention to obtaining valid area corresponding to another scheduling algorithm: Earliest Deadline (ED) (see figure 5). Referring to the previous results about the valid area obtained with a RM scheduling, the use of ED algorithm really offers new choices of task parameters : 23 points instead of 13). For example, with a strongly constraint deadline equal to 3, the period can be varied over a range as wide as 5 : P∈ [7,12]. For the choice of parameters, our server can be implemented with parameters D=3 and P=12 which unload the processor considerably while keeping the required response time for the aperiodic event.

-TEMPORAL ANALYSIS IN THE PERIOD VS RUN-TIME REPRESENTATION

In previous interpretation of the D-P plane mapping, the goal was focused on the evaluation and the improvement of the temporal characteristics of a task set. With a period versus run-time plane, the purpose is to analyze the ability of the system to undergo temporal faults (period jitters, run-time variation) without jeopardizing the schedulability of the whole task set. In the previous case, we assumed that the computational time was constant and equal to a known worst case execution time (Ci = Ci,m). But, in fact, the maximum value of the run-time of a task is quite difficult to estimate with a high accuracy because of the laborious analysis of the execution paths of the software, the compiler effects, the hardware environment (caching, paging).

Another reason to study a task with a variable runtime is that some algorithms, used in task software, contain loops that yield more accurate results with an increasing number of iterations. So, it should be very useful to know the upper limit of the execution time that a task can possess. Assuming that the deadline Di was set earlier and the possible variation range of the run-time was also determined (Ci ∈ [Ci,0, Ci,m]), we can construct the valid area of the task. First we plot the opened wedge corresponding to the initial model described in section 2.2. (zone bounded by the two lines : C=Ci,0and P=C). Then the temporal requirements of the application reduce the allowed parameters. And finally the schedulability analysis gives the final possible temporal characteristics of the task. From that result, it is quite easy to evaluate the possible variation of the execution time at a given period or the inverse. Under this framework, periodic tasks can change their execution rate or their execution times to provide different quality of service (fault tolerant, system load, etc.). This will be developed in details for a real example in the next section.

-ANALYZING A REAL EXAMPLE

-Speediness improvement of this graphical tool

A tool, based on this computer aided design methodology, was developed with a user-friendly interface that assists real-time designers in fixing the parameters of a task set of an application. And we can remark that each point of the graphical representation corresponds to a complete simulation of the application scheduling, even if the point parameters lead to a non-schedulable result. Furthermore this duration of simulation can be extremely large. So, we take care to minimize the number of parameter cases to check by using analytical conditions and properties of scheduling. We want to strongly reduce the simulated area from a number of point that is theoretically proportional to Pmax 2 in a periodic task case with a maximum period of Pmax. First, we have established some graphical properties coming from scheduling characteristics. These properties permit us to draw a line of feasible points without any schedulability calculations. For example, in RM scheduling context, we have the following results :

• property 1 : if the schedulability test is positive for the value couple (Dx, Px), it will be correct for the other value couples such as (Dx + α, Px). • property 2 : if the schedulability test is positive for the value couple (Px, Cx), it will be correct for the other value couples such as (Px, Cx -α). • property 3 : if the schedulability test is positive for the value couple (Dx, Px), it will be correct for the other value couples such as (Dx, Px + α) if the task keeps the same priority. • property 4 : if the schedulability test is positive for the value couple (Px, Cx), it will be correct for the other value couples such as (Px + α, Cx) if the task keeps the same priority.

Using these properties, associated to a dichotomic analysis among the possible couples of task parameters, the number of points to simulate is limited to number of point proportional to log(Pmax) Pmax in the RM scheduling context and to log(Pmax)+log(Pmax) with DM scheduling.

-An industrial case study

To illustrate the methodology presented in previous sections and to experiment the graphical approach, we analyze a realistic example. This industrial case study is based on an aluminum cold rolling mill application. The studied application concerns the software dedicated to the real-time check of the product quality (specially the thickness regularity). The control process software is supported and managed by another computer. The system architecture is based on both Motorola 2600 board and LynxOS™ real-time operating system. There are ten tasks presented in the Table 1 where the temporal parameters are given in milliseconds. The task set parameters used for the example denote the actual or nominal characteristics (period and execution time). The processor utilization is 0.84. The schedulability must be studied over a duration of 200 s corresponding to the LCM of task periods.

The processor is idle in the remaining time 32 s. We confine our attention here to the data processing task 3 that is constrained to be in the period range [1s, 4s], the actual period being 4 s. On the other hand, the computational time of this task can be varied from an execution to another because of the calculations depending of the data number.

Using the presented methodology, we obtain the graphical result of the schedulability domains with RM and DM scheduling algorithms (see figure 6).

For the period 2 s, the computational time can vary from 250 ms to 370 ms for a RM priority assignment and from 250 ms to 500 ms for a DM priority mapping. The major conclusion in this case is that the timing fault tolerance of the system is increased by a factor 2 concerning this task 3. If the task period can be set to the upper value of the specified period range (P3,max = 4 s), the execution time is in the range [250 ms, 750 ms] for RM scheduling and [250 ms, 1.1 s] for DM scheduling.

In this latter case (P3,max = 4 s, DM scheduling), the application is schedulable for a processor utilization varying from 0.78 to close to 1 due to the task 3 execution time variation. Fig. 6 : The allowed temporal parameters of the task 3 of the example described in Table 1.

-CONCLUSION

In this paper, we have introduced a graphical representation for tasks based on a classical temporal model. The elaboration of the valid zone for a given task is computed in three steps which take into account the intrinsic characteristics of the periodic task, the timing requirements of the application for this particular task and finally the schedulability of the task set. These successive surfaces are of course smaller and smaller. At the end, the user has to take design decisions and to chose the right parameters of the task in the way to improve the temporal performances of the task itself and more globally of the application.

With a representation in deadline-period co-ordinate plane, this approach is particularly useful in situations in which the designer has to determine parameters of a task set, either to add a task to a set of tasks with already fixed parameters or to modify the task parameters in order to analyze and improve some temporal characteristics (response time, processor load, …). With a specific representation (period versus run-time plane), the purpose is to analyze the ability of the system to tolerate temporal faults (period jitters, run-time variation, …) without risking the schedulability of the whole task set. Moreover, we have shown how this graphical model permits us to compare the different scheduling algorithms : RM, DM and ED.

Fig. 1 :

 1 Fig. 1 : A graphical representation of the possible values of temporal parameters of a task.

Fig. 2 :

 2 Fig. 2 : A graphical representation of the possible values of temporal parameters (Di, Pi) for periodic tasks corresponding to the temporal requirements.

Fig. 3 :

 3 Fig.3: The allowed temporal attributes of the periodic server in a multitasking environment.

Fig. 5 :

 5 Fig. 4 : Example with a RM scheduling.

Table 1 :

 1 A real-time application with ten tasks

	Task	C (ms)	D (s)	P (s)
	T1	250	1	2
	T2	500	4	4
	T3	250	4	4
	T4	500	4	4
	T5	250	4	4
	T6	500	3	8
	T7	1,250	10	20
	T8	2,500	50	100
	T9	12,500	200	200
	T10	12,500	200	200


Camille Le Roy
email: camille.leroy@wur.nl

Dario Amadori

Samuel Charberet

Jaap Windt

Florian Muijres

Violaine Llaurens

Vincent Debat

Adaptive evolution of flight in Morpho butterflies

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

Abstract:

The diversity of flying animals suggests that countless combinations of flight morphologies and behaviors have evolved with specific lifestyles, thereby exploiting diverse aerodynamic mechanisms. Elucidating how morphology, flight behavior and aerodynamic properties together diversify with contrasted ecology remains rarely accomplished. Here, we studied the adaptive co-divergence in wing shape, flight behavior and aerodynamic efficiency among Morpho butterflies living in different forest strata, by combining high-speed videography in the field with morphometric analyses and aerodynamic modelling. By comparing canopy and understory species, we show that adaptation to an open canopy environment resulted in increased glide efficiency. Moreover, this enhanced glide efficiency was achieved by different canopy species through strikingly distinct combinations of flight behavior, wing shape and aerodynamic mechanisms, highlighting the multiple pathways of adaptive evolution.

One Sentence Summary: Amazonian Morpho butterflies reveal a combined adaptive divergence of wing shape and flight behavior across forest strata.

Main Text:

Insects display a diversity of flight patterns reflecting their different ecologies, from the sustained, energy-efficient flight of long-range migrating species [START_REF] Gibo | Insect Flight[END_REF], to the highly-maneuverable hovering Sympatric Morpho species display remarkably contrasted ecologies, some species flying in the dense vegetation of the understory and others in the open canopy [START_REF] Devries | Vertical distribution, flight behaviour and evolution of wing morphology in Morpho butterflies[END_REF]. The divergence between canopy and understory species occurred around 22 million years ago [START_REF]Punctual ecological changes rather than global factors drive species diversification and the evolution of wing phenotypes in Morpho butterflies | bioRxiv[END_REF]. This resulted in contrasted selective pressures acting on the evolution of flight behavior and morphology among species, whereby open canopy habitats may favor a more extensive gliding behavior [START_REF] Devries | Vertical distribution, flight behaviour and evolution of wing morphology in Morpho butterflies[END_REF]. We examined the divergence of flight behavior among habitats and asked whether the evolution of gliding flight in canopy species was enabled by increased aerodynamic efficiency through changes in wing shape [START_REF] Chazot | Morpho morphometrics: Shared ancestry and selection drive the evolution of wing size and shape in Morpho butterflies[END_REF].

We performed a series of field and semi-field experiments in Amazonian Peru. Here, up-totwelve Morpho species co-occur, allowing us to investigate how habitat impacts the evolution of flight in closely-related species in sympatry. We used high-speed videography to track and characterize flight behavior of wild individuals in the field and in a large insectary (Movies S1-S6). We then quantified the wing shape of these filmed butterflies using geometric morphometrics and assessed shape covariation with flight. Finally, we used computational fluid dynamics (CFD) modelling to assess the aerodynamic efficiency associated with the contrasted wing shapes of species specialized in different habitats.

Using high-speed videography, we recorded 136 sequences of 80 wild Morpho butterflies freely patrolling in nature, including four understory and three canopy species (Fig. 1, Movies S3-S6).

From the temporal positions of each wing stroke, we measured the flapping frequency, the gliding phase duration, and the temporal flap-gliding ratio of each flight (Fig. 1). While little variation was found in flapping frequency, two of the canopy species (Morpho cisseis and Morpho telemachus) showed sharply longer gliding phases than all understory species, spending about half of their time gliding. The third canopy species, Morpho rhetenor, differed from the other canopy species, showing limited use of gliding, even less than all understory species (Fig. 1). These findings corroborate field observations [START_REF] Devries | Vertical distribution, flight behaviour and evolution of wing morphology in Morpho butterflies[END_REF] and highlight that contrasted flight behaviors exist within the canopy clade.

To finely characterize flight behaviors of canopy and understory species, we built a large outdoor insectary equipped with a high-speed stereoscopic videography system. In this large cage, we tracked the three-dimensional movements of 241 flights from 82 wild-caught Morpho butterflies from eight understory and three canopy species (Fig. 2A,B, Movies S1-S2). We then characterized all flights using eleven flight kinematics parameters, of which six characterized the complete flap-gliding flight, three the flapping-phase, and two the gliding phase (supplementary text, Table S1).

We first compared insectary flights to those recorded in the wild to assess the impact of captivity on flight behavior (Fig. 1, Table S2). Captivity reduced the extent of flap-gliding flight in canopy species. Understory species were little affected, suggesting that they may be accustomed to flying in confined spaces. Overall, the interspecific variation in flight behavior in the insectary was broadly consistent with that observed in the wild (Fig. 1).

A Principal Component Analysis (PCA) performed on the eleven kinematics parameters showed that flight behaviors significantly differed between canopy and understory species (Fig. 2C). A phylogenetic MANOVA confirmed that this difference is higher than expected from a Brownian model of character evolution (Table S3). This strong divergence in flight mode between canopy and understory species therefore cannot be explained by their phylogenetic divergence only, therefore pointing at an effect of the contrasted selection regimes acting on flight evolution in the two microhabitats.

Principal component 1 (PC1) was driven by the relative use of gliding flight (variation in glideduration, glide-angle and glide-ratio), which was comparable in both canopy and understory species when flying in captivity (Fig. 2C). PC2 reflected the aerodynamic force production during flapping flight, and clearly opposed canopy and understory species: fast flight and high advance ratio for understory species on the negative values and slow flight and curvy trajectories for canopy species on the positive values. Understory species thus exhibit a more powerful wingbeat, producing higher aerodynamic forces and leading to higher advance ratios and straighter highspeed flights.

Glide angle and ascent angle diverged more between canopy and understory species than predicted phylogenetically (Fig. 2D, Table S3), suggesting a strong effect of natural selection on these two flight components. During their few flapping phases, canopy butterflies also climbed more steeply, their mean ascent angle being 70% larger than in understory species (Fig. 2D; ascend,canopy=22°±3°, mean±standard deviation, n=70 flapping phases; ascend,understory=13°±1°, n=171 flapping phases). Thus, although understory species tended to fly at higher advance ratios and flight speeds (Fig. 2C,PC2), the ascent angle of canopy species was higher (Fig. 2C,PC1; Fig. 2D). This could stem from an increased behavioral tendency to fly up, and/or from a higher climbing efficiency generated by their morphology.

During the gliding phases, canopy butterflies had a 36% lower glide angle than the understory butterflies (Fig. 2D; glide,canopy=7°±2°, n=61 gliding phases; glide,understory=11°±1°, n=135 gliding phases). Such shallow glides allow canopy butterflies to travel longer distances for a given height loss, consistent with the longer gliding phases measured in the wild. Glide angle is directly related to the aerodynamic efficiency parameter lift-to-drag ratio [START_REF] Norberg | Structure, form, and function of flight in engineering and the living world[END_REF], and shallower angles detected in canopy species might be promoted by their divergent wing shapes. This combination of field and semi-field experiments shows that the evolutionary shift from understory to canopy resulted in an increased use and efficiency of gliding flight (Fig. 2D), combined with a reduction in aerodynamic force production during forward flapping flight (Fig. 1C,PC2).

We then investigated the contribution of morphological divergence in the adaptive evolution of flight between canopy and understory species. We precisely quantified wing shape of the filmed butterflies using geometric morphometrics and detected a strong covariation between wing shape and flight behavior using a phylogenetic partial least square analysis (Fig. 3A). Butterflies with more rounded wings and higher wing-loading flew at higher flight speed and advance-ratio, and accelerated more rapidly. These results suggest that the evolution of smaller (high wing-loading), more-rounded wings indeed increased force production during flapping flight. Our analyses demonstrate that flight power therefore tightly co-evolves with wing shape.

In contrast to flapping flight parameters, gliding parameters were weakly correlated with wing shape and aspect-ratio AR (Fig. 3B, Table S4). Canopy species are efficient gliders yet exhibit strikingly diverse wing-loadings, aspect-ratios and wing shapes (Fig. 3): two of the studied canopy species are slow flyers with low aspect-ratio, low wing-loading triangular wings (M. cisseis and M. theseus), whereas the fast-flying M. rhetenor has high aspect-ratio elongated wings with high wing-loading. This begs the question of how this divergence in wing shape among species altered gliding efficiency.

Using Computational Fluid Dynamics (CFD), we then determine how glide performance differs between canopy and understory species (ncanopy=4; nunderstory=3; Figs. 1,4). For each species, we produced in-silico wings based on our gliding flight experiments (Figs. 4E-J,S2). We then performed gliding-flight CFD simulations to determine the lift-to-drag-ratio to angle-of-attack curves (L/D- Fig. 4A). Maximum lift-to-drag-ratio L/Dmax was achieved at °°, and was 9% greater in canopy species (L/Dmax,canopy=5.62±0.19; L/Dmax,understory=5.18±0.07; phylogeneticgeneralized-least-squares: F1,5=13.48, P=0.01; Fig. 4A,B), indicating that their wing shapes confer higher glide efficiency.

Wing shape primarily affects induced drag of a wing, which inversely scales with the product of wing aspect-ratio AR and span-efficiency e [START_REF] Muijres | Actuator disk model and span efficiency of flapping flight in bats based on time-resolved PIV measurements[END_REF]. Therefore, we tested how L/Dmax scaled with these parameters, and how this varied between species (Fig. 4BD). The product AR•e was 24% higher for canopy species (AR•ecanopy=1.83±0.15; AR•eunderstory=1.48±0.05; phylogeneticgeneralized-least-squares: F1,5=15.38, P=0.01; Fig. 4B), which was achieved by canopy species in different ways (Figs. 4C,D): M. rhetenor has exceptionally high aspect-ratio wings (Figs. 4C), whereas the other species have primarily an enhanced span-efficiency (Fig. 4D). Airflow visualizations at L/Dmax (Fig. 4E-J) show that all gliding butterflies produce a stable Leading-Edge-Vortex and streamwise wingtip-vortices. Unlike canopy species, all understory species also produce a highly-turbulent wingroot-vortex, which explains their reduced aerodynamic efficiency [START_REF] Muijres | Actuator disk model and span efficiency of flapping flight in bats based on time-resolved PIV measurements[END_REF]. These results provide functional evidence that wing shape divergence among Morpho species directly affects glide efficiency, and that different canopy species with contrasted wing shapes achieve this using distinct aerodynamic properties. Divergence in wing shape across forest strata has been documented in Amazonian butterflies [START_REF] Devries | Vertical distribution, flight behaviour and evolution of wing morphology in Morpho butterflies[END_REF][START_REF] Chazot | Morpho morphometrics: Shared ancestry and selection drive the evolution of wing size and shape in Morpho butterflies[END_REF][START_REF] Cespedes | Cruising the rain forest floor: butterfly wing shape evolution and gliding in ground effect[END_REF]. Because any trait that reduces energetic cost-of-flight is likely under positive selection [START_REF] Norberg | Structure, form, and function of flight in engineering and the living world[END_REF], evolution in open habitats such as the canopy may favor traits enhancing glide efficiency. Most animals flying in open environments indeed display this energy-saving gliding behavior [START_REF] Gibo | Insect Flight[END_REF][START_REF] Rayner | Current ornithology[END_REF].

Our combination of aerodynamic and ecological approaches revealed

Intriguingly, our study also reveals an unexpected flight behavior in a canopy species that mostly uses flapping flight (M. rhetenor; Fig. 1 and Movies S4). This discrepancy between aerodynamic performance and behavior suggests that conflicting selective pressures affect flap-gliding behavior. The vigorous flight of M. rhetenor might have co-evolved with its blue iridescence, as the blue flashes induced during wing flapping cause confusion in predators [START_REF] Murali | Now you see me, now you don't: Dynamic flash coloration as an antipredator strategy in motion[END_REF]. Despite these differences in flight behavior, all canopy species show increased glide efficiency compared to understory species (Fig. 4), suggesting that the selection of aerodynamically-efficient wing shapes prevails. Overall, our study illustrates how adaptive evolution is fueled by the flexible adjustment of morphology, behavior and aerodynamic performance.

  of nectar-feeding species (2). This diversity of flight modes stems from morphological and behavioral adaptations that improve flight performance metrics such as speed, maneuverability or energetic efficiency. Investigating insect flight aerodynamics is therefore crucial to understand how natural selection shapes flight. Although insect flight has been studied in detail in several species, including Drosophila, mosquitoes and hawkmoths (3), only the comparison of closelyrelated species adapted to different habitats can unravel the impact of ecological constraints on the diversification of aerodynamic properties. Here, we addressed the ecological, behavioral and morphological bases of the diversification of flap-gliding flight in closely-related butterfly species. Butterflies are the only insects that regularly use flap-gliding flight, by combining periods of flapping interspersed with gliding. In contrast, many intermediate-sized bird species use flap-gliding flight to reduce energetic expenditure when the aerodynamic efficiency of gliding phases is high enough (4). In this study, we assessed the diversity of flap-gliding flight in the neotropical butterfly genus Morpho.

  how natural selection imposed by different microhabitats can drive the evolution of flap-gliding flight by jointly altering wing shape and flight behavior. Butterflies from species evolving in the cluttered understory habitat display powerful flapping phases, resulting in high flight speeds and advance ratios. In contrast, evolution in the open canopy resulted in a more efficient gliding flight, illustrated by the reduced descend angles during gliding phases observed in the canopy species.

Fig. 1 .

 1 Fig. 1. Canopy butterflies use flap-gliding flight in a larger extent than understory butterflies. The twelve studied Morpho species are represented with their phylogenetic relationships. Differences in flap-gliding parameters between microhabitats were more striking in nature, as captivity reduced gliding in canopy species. Bars indicate the mean ± standard error, and stars indicate significant difference between nature and captivity. Note that M. helenor and M. achilles cannot be distinguished during flight: corresponding data in nature were thus pooled.

Fig. 2 .

 2 Fig. 2. Flight kinematics revealed differences in behavior and performance between canopy and understory species. (A) The 241 trajectories analyzed are shown together. (B) A single flight trajectory of an individual M. cisseis (duration=1.7 sec). Droplets indicate the uppermost and lowermost wing positions during up-stroke and downstroke, respectively. (C) Principal Component Analysis showing the divergence of flight between canopy and understory butterflies.(D) Gliding and climbing efficiency was higher in canopy species, and found to diverge more strongly than expected from phylogenetic distance.

Fig. 3 .

 3 Fig. 3. Wing shape and wing-loading jointly covary with flight behavior. (A) Phylogenetic partial-least square analysis shows the covariation between wing shape and flight behavior (r-PLS = 0.89; P = 0.02; 74% of covariation explained). This covariation opposes triangular to rounded wings, respectively associated with slowcurvy flight, and straightermore powerful flight. Wing-loading (depicted by symbol size) also covaries with flight, suggesting that the evolution of flight is linked to both wing shape and body morphology. Flight loadings are indicated on the right. (B) Phylogenetic morphospace depicting variation in wing shape among species, diamond size indicates aspect ratio.

Fig. 4 .

 4 Fig. 4. Canopy species achieve greater glide efficiency through different wing shapes. (A,B) Canopy species have increased maximum lift-to-drag ratio (L/Dmax) and aerodynamic efficiency (AR•e), as shown by phylogenetic-generalized-least-squares analyses (top and right of B, respectively). (B-D) Canopy-species M. rhetenor and M. cisseis achieve increased aerodynamic efficiency via increased aspect ratio (C) and span efficiency (D), respectively. (E-J) Morpho wings at L/Dmax, with relative air pressure and streamlines (left wing) and vorticity fields color-coded with turbulence state (right wing). Understory species produce highly-turbulent wingroot vortices (E-G), canopy species do not (H-J). Species-specific error bars show numerical uncertainties, and for canopy/understory groups these are standard errors.

Acknowledgments:

We thank the Peruvian authorities, and in particular SERFOR for providing the necessary research permits (permit: 002-2015-SERFOR-DGGSPFFS). We thank Ronald Mori Pezo and Joel Pintado for help with butterfly sampling. We thank Ty Hedrick and Guilherme Vaz for their advice on the videography and CFD analyses, respectively. Funding: This work was supported by grants from the Agence National de la Recherche under the LabEx ANR-10-LABX-0003-BCDiv, in the program "Investissements d'avenir" ANR-11-IDEX-0004-02 (to C.L.R), the Emergence program of Paris city council (to V.L.), and Université de Paris and the Ecole Doctorale

FIRE -Program Bettencourt (to C.L.R.). Author contributions: C.L.R., V.L., and V.D. designed the research plan and performed the field experiments. S.C. analyzed experimental field data. D.A., F.T.M., and J.W. designed and performed the CFD study. C.L.R., V.L., V.D. and F.T.M. wrote the paper. All authors edited the manuscript. Data and materials availability: All data relevant to the manuscript are provided in the supplementary materials.

Supplementary Materials:

Materials and Methods

Figs. S1 to S13

Tables S1 to S6

Captions for Movies S1 to S6 Source data (excel file)