Mikel Cordovilla

Frédéric Boniol

Julien Forget

Eric Noulard

Claire Pagetti

Developing critical embedded systems on multicore architectures: the PRELUDE-SCHEDMCORE toolset

In this paper we present an end-to-end framework for the design and the implementation of embedded systems on a symmetric multicore. The developer first specifies the system using the PRELUDE language, a formal realtime architecture description language. The compiler then translates the program into a set of communicating periodic tasks that preserves the semantics of the original program. The schedulability analysis is performed by the SCHEDMCORE analyzer. If the program is schedulable, it can finally be executed on the target multicore architecture using the SCHEDMCORE execution environment.

Introduction

The development of critical embedded systems undergoes strict development processes. In this context, the purpose of this paper is to propose an integrated development framework for critical embedded systems, that goes from the programming of high level specifications, using the PRELUDE language, to the execution of a semantically equivalent multithreaded code (generated automatically) on a symmetric multicore architecture, using the SCHEDMCORE framework.

Targeted systems

As an example, let us consider an adapted version of the Flight Application Software (FAS) of the Automated Transfer Vehicle (ATV) designed by Astrium for resupplying the International Space Station (ISS). Figure 1 provides a simplified informal description of its software architecture. The FAS acquires several data treated by dedicated sub-functions: orientation and speed (Gyro Acq), position (GPS Acq and Str Acq) and telecommands from the ground station (TM/TC). The Guidance Navigation and Control function (divided into GNC_US and GNC_DS) computes the commands to apply while the Failure Detection Isolation and Recovery function (FDIR) verifies the state of the FAS and checks for possible failures. Commands are sent to the control devices: thruster orders (PDE), power distribution orders (PWS), solar panel positioning orders (SGS) and telemetry towards the ground station (TM/TC).

Figure 1. FAS architecture

The specification of such a system must precisely define three things: [START_REF] Aubry | Synchronous distribution of Signal programs[END_REF] the functional behaviour of each function: each function can be programmed separately with existing languages, such as SIMULINK, SCADE or directly in low level code (C for instance); (2) the real-time characteristics of the system: periods and deadlines which depend on the physical characteristics of the system; (3) multi-rate communication patterns: the way data is exchanged between functions is complex since functions of different periods communicate.

A correct implementation must respect all the real-time constraints of the system and must also be functionally deterministic, meaning that the outputs of the system must always be the same for a given sequence of inputs. This requires task communications to be fully deterministic, which means that for any execution of the system the same job (or instance) of the producer must communicate with the same job of the consumer of the communication.

Contribution

Figure 2 shows the development process within the PRELUDE-SCHEDMCORE framework. In this paper, we focus more on presenting a complete operational toolset rather than on new theoretical results (we published most results separately previously). We show how the different steps of the development come together inside our framework and we detail the tools developed to this intent. The framework is made up of three components:

• The PRELUDE compiler which has been extended to generate code executable by the SCHEDMCORE environment (Section 3);

• The SCHEDMCORE multiprocessor schedulability analyzer (Section 4);

• The SCHEDMCORE environment which is an extensible, easy-to-use and portable real-time scheduler framework for multicore (Section 5).

PRELUDE specification

i m p o r t e d node t a u k ( i : i n t ) r e t u r n s ( j : i n t ) w c e t 7 ; node c o n s o ( i ) r e t u r n s ( o ) l e t o= t a u 2 ( ( 0 f b y t a u 1 ( i ) ) * ˆ4 ) ; t e l SCHEDMCORE libraries FP, gEDF, LLREF, off line PRELUDE [START_REF] Forget | A multi-periodic synchronous data-flow language[END_REF][START_REF] Forget | A Synchronous Language for Critical Embedded Systems with Multiple Real-Time Constraints[END_REF] 1 is a formal language designed for the specification of the software architecture of a critical embedded control system. It belongs to the synchronous data-flow languages [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] and focuses on dealing with the functional and real-time aspects of multi-periodic systems conjointly. From a PRELUDE program the compiler generates a program consisting of a periodic dependent task set. For each pair of producer job and consumer job of the task set the following properties must be ensured: (1) the producer completes before the consumer starts. This is modeled by adding a precedence constraint from the producer to the consumer. The respect of the precedence constraints must be ensured by the scheduler. (2) Produced data remains available until the completion of the consumer. This is fulfilled using a specific communication protocol (directly generated by the compiler) derived from [START_REF] Sofronis | A memory-optimal buffering protocol for preservation of synchronous semantics under preemptive scheduling[END_REF] and detailed in [START_REF] Forget | A Synchronous Language for Critical Embedded Systems with Multiple Real-Time Constraints[END_REF]. The initial release of the compiler (described in [START_REF] Forget | A Synchronous Language for Critical Embedded Systems with Multiple Real-Time Constraints[END_REF]) was targeted for a uniprocessor platform. The task set was executed with MARTE OS [START_REF] Rivas | A POSIX-Ada Interface for Application-Defined Scheduling[END_REF], using a scheduling policy derived from EDF [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard-real-time environmet[END_REF] and from the work of [START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF]. In this paper we present how it was extended to enable the execution of PRELUDE programs on a multicore architecture.

Multiprocessor schedulability analyzer

The task set generated by PRELUDE must be scheduled in a way that respects the real-time attributes of each task and the extended precedence constraints between the tasks (the constraints can relate tasks with different periods).

Scheduling policies on multiprocessor platform have been summarized in surveys such as [START_REF] Davis | A survey of hard real-time scheduling algorithms and schedulability analysis techniques for multiprocessor systems[END_REF]. We consider preemptive policies accepting (full) migration, meaning that a suspended instance may be resumed on a different processor. The SCHEDMCORE2 framework implements the policies FP (fixed priority), gEDF (global Earliest Deadline First), gLLF (global Least Laxity First) and LLREF (Largest Local Remaining Execution First) [START_REF] Cho | An optimal real-time scheduling algorithm for multiprocessors[END_REF].

Before executing the program, we must check that the execution will respect the real-time constraints of the system. This schedulability analysis determines whether a task set is correctly scheduled by a given policy on an architecture. We have developed an exact approach based on an efficient enumerative technique [START_REF] Cordovilla | Multiprocessor schedulability analyser[END_REF] for dependent task sets. In summary, the possible execution sequences of a system are first encoded as a finite automaton. Then, using either the UPPAAL model checker [START_REF] Behrmann | A tutorial on UPPAAL[END_REF] or a custom C program, we explore this automaton to find which executions violate the constraints of the system. Though this brute force method may seem expansive, we obtain pretty good results even with complex task sets (huge number of tasks, utilization factor close to the number of processors or huge least common multiple of periods). The automaton model can also be used to generate a schedule computed off-line.

SCHEDMCORE execution environment

The objective of SCHEDMCORE execution environment (or runtime) is to provide a validation and execution tool, somewhere between a pure simulator and a true hard realtime execution environment. It allows to run a set of tasks written in C using various real-time multicore scheduling policies on a standard multicore system. Most simulators are dedicated either to functional simulation or temporal simulation, while SCHEDMCORE enables to do both simultaneously. Even if SCHEDMCORE runtime cannot be considered as a hard real-time runtime (yet), it is far more realistic than a simulator: it can schedule real C functions on a real-time timescale. This provides a "simulation" pretty faithfull to the actual system, while remaining portable, easy to use and to extend. SCHEDMCORE provides libraries for most common multicore scheduling policies but was designed so that a user can easily add his own new policy.

Related works

Development frameworks for multi-periodic systems SIMULINK [START_REF]Simulink: User's Guide[END_REF] is widely used in many industrial appli-cation domains and allows the description of communicating multi-periodic systems. Current code generators, such as Real-time workshop-embedded coder from the MathWorks or TargetLink from dSpace, provide a multithreaded code generation on uniprocessor. In [START_REF] Caspi | From simulink to scade/lustre to tta: a layered approach for distributed embedded applications[END_REF] the authors translate a multi-periodic SIMULINK specification into a semantical equivalent synchronous program that is executed on a multiprocessor time triggered architecture (TTA) where no preemption is allowed.

Synchronous data-flow languages, such as LUSTRE [START_REF] Halbwachs | The synchronous data-flow programming language LUSTRE[END_REF] or SCADE [START_REF] Dormoy | Scade 6 a model based solution for safety critical software development[END_REF], have been extended with operators that enable the specification of real-time constraints in order to program multi-threaded system more easily [START_REF] Curic | Implementing Lustre programs on distributed platforms with real-time constraints[END_REF]. Thread synchronizations rely on a specific communication protocol initially defined in [START_REF] Sofronis | A memory-optimal buffering protocol for preservation of synchronous semantics under preemptive scheduling[END_REF] for uniprocessor and later extended specifically for Loosely Time Triggered Architectures (LTTA) in [START_REF] Tripakis | Implementing synchronous models on loosely time-triggered architectures[END_REF].

Finally, automated distribution of synchronous program has been studied in [START_REF] Girault | Automatic rate desynchronization of embedded reactive programs[END_REF][START_REF] Aubry | Synchronous distribution of Signal programs[END_REF]. However, these studies are not dedicated to multi-periodic systems and thus scheduling policies are not considered.

Schedulability analysis Lots of theoretical results are already available for independent task sets (see surveys [START_REF] Baker | Schedulability analysis of multiprocessor sporadic task systems[END_REF][START_REF] Davis | A survey of hard real-time scheduling algorithms and schedulability analysis techniques for multiprocessor systems[END_REF]). However there are not so many for dependent task sets and not so many tools (even for independent task sets) are available yet. STORM [START_REF] Urunuela | STORM, simulation tool for real-time multiprocessor scheduling[END_REF] is a multiprocessor scheduling simulator. In [START_REF] David | Model-Based Design for Embedded Systems, chapter Model-Based Framework for Schedulability Analysis Using UP-PAAL 4[END_REF] a framework allowing the analysis of tasks configuration for multiprocessor machines with UPPAAL models is proposed. This framework supports rich task models including timing uncertainties in arrival and execution times, and dependencies between tasks (such as precedences and resource occupation). However the task set should not exceed 30 tasks due to performance concerns.

Execution environments Lots of real-time execution environments or operating systems are available: industrial ones like VxWorks, LynxOS or PikeOS; modified Linux like Xenomai, LitmusRT or SCHED EDF for Linux (see [START_REF] Faggioli | An edf scheduling class for the linux[END_REF] and references therein); academics OS like MarteOS [START_REF] Rivas | A POSIX-Ada Interface for Application-Defined Scheduling[END_REF] and even user space runtime like Metascheduler [START_REF] Li | A formally verified application-level framework for real-time scheduling on posix real-time operating systems[END_REF]. We need an open environment that enables the implementation of user-specific scheduling policies, which makes us rule out industrial solutions. Modified Linux like LitmusRT or sub-kernel approaches like Xenomai are too closely tied to the kernel for our purpose as this requires to patch and recompile the kernel when new releases of Linux occur. SCHEDMCORE can run on top of these environments but does not rely on them. MARTE OS [START_REF] Rivas | A POSIX-Ada Interface for Application-Defined Scheduling[END_REF] and the Meta-scheduler [START_REF] Li | A formally verified application-level framework for real-time scheduling on posix real-time operating systems[END_REF] answer our needs partially but, in both cases, those environments only support uniprocessor platform. In the end, we reused the conceptual idea of a pluggable scheduler framework (of MARTE OS and the Meta-scheduler) but started with a fresh new source code.

The PRELUDE language

For the first step of the development, that is the specification of the system, we rely on the PRELUDE language previously introduced in [START_REF] Forget | A multi-periodic synchronous data-flow language[END_REF][START_REF] Forget | A Synchronous Language for Critical Embedded Systems with Multiple Real-Time Constraints[END_REF]. PRELUDE is a formal Real-Time Architecture Description Language with synchronous semantics dedicated to the specification of critical embedded control systems.

Language definition

Flows and clocks PRELUDE is a synchronous data-flow language and thus shares similarities with LUSTRE [START_REF] Halbwachs | The synchronous data-flow programming language LUSTRE[END_REF], SIGNAL [START_REF] Benveniste | Synchronous programming with events and relations: the Sig-nal language and its semantics[END_REF] or LUCID SYNCHRONE [START_REF] Pouzet | Lucid Synchrone, version 3. Tutorial and reference manual[END_REF]. Variables and expressions are flows. A flow is a sequence of pairs (v i , t i ) i∈N , where v i is a value in some domain V and t i is a date in Q (∀i, t i < t i+1 ). The clock of a flow, is its projection on Q. It defines the set of instants during which the values of the flow are computed: value v i must be computed during the interval (or instant) [t i , t i+1 [. According to this relaxed synchronous hypothesis, different clocks have different durations for their instants. Two flows are synchronous if they have the same clock. We distinguish a specific class of clocks corresponding to periodic task activations called strictly periodic clocks: the clock h = (t i ) i∈N is strictly periodic if there exists some n such that t i+1 -t i = n for all i. n is the period of h and t 0 is its phase.

Nodes hierarchy A PRELUDE program consists of a set of nodes. A node defines its outputs flows from its input flows. It can either be hierarchical, in which case the relationship between its inputs and outputs is defined by a set of equations, or it can be imported, in which case the node is implemented outside the program using another existing language (C or LUSTRE for instance). A simple program is given below: imported node p l u s 1 ( i : i n t ) r e t u r n s ( o : i n t ) wcet 5 ; node N( i : i n t r a t e ( 1 0 , 0 ) ) r e t u r n s ( o : i n t r a t e ( 1 0 , 0 ) ) l e t o= p l u s 1 ( i ) ; t e l Following the data-flow programing style, the node N is activated each time it receives an input i, i.e. each 10 time units (the term (10,0) denotes the clock of period 10 and of phase 0). At each activation, it computes the output o as the result of the application of the imported node plus1 to i.

Rate transition operators PRELUDE defines a set of rate transition operators that enable the definition of userspecified communication patterns between nodes of different rates. Let e be a flow expression, k ∈ N and q ∈ Q:

• cst fby e is the delay operator borrowed from LU-CID SYNCHRONE: it first produces cst and then produces each value of e delayed by one instant. The clock of this flow is the same as that of e;

• e * ˆk produces a flow whose period is k times shorter than that of e. Each value of e is duplicated k times in the result;

• e/ ˆk produces a flow k times slower than e. The result only keeps the first value of each k successive values of e;

• e ∼> q produces a flow where each value of e is delayed by q * n, where n is the period of e.

These operators are illustrated on the example below (nodes τ 1 , τ 2 and τ 3 are imported, their declaration is omitted here):

node s a m p l i n g ( i : r a t e ( 1 0 , 0 ) ) r e t u r n s ( o1 , o2 )

var vf , v s ; l e t ( o1 , v f ) = t a u 1 ( i , ( 0 f b y v s ) * ˆ3 ) ; v s = t a u 2 ( v f / ˆ3 ) ; o2 = t a u 3 ( ( v f ˜> 1 / 1 0 ) / ˆ6 , ( v s ˜> 1 / 3 0 ) / ˆ2 ) ; t e l
The behaviour of this program is the following (we detail the values produced and their dates of production): date 0 1 10 20 30 40 50 60 61 ...

vf v 0 v 1 v 2 v 3 v 4 v 5 v 6 ... vf/ˆ3 v 0 v 3 v 6 ... vs v ′ 0 v ′ 1 v ′ 2 ... 0 fby vs 0 v ′ 0 v ′ 1 ... (0 fby vs) * ˆ3 0 0 0 v ′ 0 v ′ 0 v ′ 0 v ′ 1 ... vf∼>1/10/ˆ6 v 0 v 6 ... vs∼>1/30/ˆ2 v ′ 0 v ′ 2 ...

Compilation

The compiler translates a PRELUDE program into a set of dependent periodic tasks. The preservation of the semantics of the program through the different steps of the compilation was formally proved in [START_REF] Forget | A Synchronous Language for Critical Embedded Systems with Multiple Real-Time Constraints[END_REF].

Static analyses

The compiler first performs a series of static analysis, to ensure that the program is correct, i.e. that its semantics is well-defined. First a standard type-checking is performed [START_REF] Pierce | Types and programming languages[END_REF]. The language is strongly typed, in the sense that the execution of a program cannot produce a run-time type error. Then, the causality check (derived from [START_REF] Halbwachs | Generating efficient code from data-flow programs[END_REF]) verifies that the program is causal, which means that no variable instantaneously depends on itself (not without a fby in the dependencies). Finally, the clock calculus [START_REF] Forget | A multi-periodic synchronous data-flow language[END_REF] verifies that the program only combines flows that have the same clock. This ensures that the program only accesses values at dates at which they are indeed available.

Translation into a dependent task set

The system model produced by PRELUDE consists of a set of concurrent periodic tasks communicating via shared variables stored in buffers. This model is encoded in a C file generated by the compiler. More formally, a system is defined as a tuple S, V, R, C where:

• S = {τ i } i=1,...,n is a finite periodic task set;
• V is the set of variables produced and consumed by the tasks. τ i .in ⊆ V (resp. τ i .out ⊆ V) stands for the set of variables that are consumed (resp. produced) by τ i ;

• R ⊆ S ×P ×S is the precedence relation, where P is the set of periodic extended precedence constraints; ). This task model is a standard adaptation from [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard-real-time environmet[END_REF].

• C : S × S × V → W
Extended precedence constraints Tasks are dependent because they are related by precedence constraints. A simple precedence constraint τ i → τ j relates two tasks with the same period and imposes that each job of the producer executes before one job of the consumer. An extended precedence constraint relates a subset of the jobs of the communicating tasks. Let τ i.n → τ j.n ′ denote a precedence constraint from the instance n of τ i to the instance n ′ of τ j . For any n ∈ N, let I n denote the set of integers of the interval [0, n[.

Definition 1 (Periodic Extended Precedence Constraint)

Let τ i and τ j be two tasks, p = lcm(T i , T j ) and M i,j ⊆ I p/Ti × I p/Tj (notice that M i,j is a finite set).

The periodic extended precedence constraint τ i Mi,j ---→ τ j is defined as the following set of task instance precedence constraints:

∀(n, n ′ ) ∈ M ′ i,j , τ i.n → τ j.n ′ M ′ i,j = {(n, n ′ )| ∃k ∈ N, (m, m ′ ) ∈ M i,j , (n, n ′ ) = (m, m ′ ) + (k p Ti , k p Tj ) }
A simple precedence constraint τ i → τ j is actually a particular case of periodic extended precedence constraints where M i,j = {(0, 0)}. Other examples are drawn below.

Communication protocol Tasks are synchronous meaning that all their inputs τ i .in must be available at the beginning of the execution of each job and all their outputs τ i .out are produced simultaneously at the end of the execution of each job.

There is a data-dependency between two tasks τ i and τ j when there exists a variable v that is produced by τ i and τi τi τi τi τi τi τj τj

(a) M i,j = {(0, 0)}

τ i τ i τ j τ j τ j τ j τ j τ j (b) M i,j = {(0, 2)} τi τi τi τi τi τj τj τj (c) M i,j = {(0, 0), (2, 1), (4, 3)} τ i τ i τ j τ j τ j (d) Not representable
consumed by τ j . The compiler allocates a finite buffer b i,j,v where τ i writes values and τ j reads values. The protocol specifying how a task accesses the buffer is defined by a dependence word. Let the ultimately periodic word w = v(u) ω denote the infinite sequence of integers consisting of the finite prefix v followed by the infinite repetition of the finite sequence u. Let w[k] denote the k th value of w.

Definition 2 (Dependence word) Let τ i , τ j be two tasks and v be a variable such that v ∈ τ i .out and v ∈ τ j .in. The dependence words w = C(τ i , v, τ j ) and w ′ = C(τ j , v, τ i ) are ultimately periodic words such that w[k] denotes the index of the cell of b i,j,v where τ i.k writes and w ′ [k] denotes the cell where τ j.k reads.

For instance, if C(τ i , v, τ j ) = 0(102) ω , then values produced by τ i are stored in the buffer as follows: τ i.0 is not stored, τ i.1 is stored in cell 1, τ i.2 is not stored, τ i.3 is stored in cell 2, then τ i.4 is stored in cell 1, τ i. [START_REF] Benveniste | Synchronous programming with events and relations: the Sig-nal language and its semantics[END_REF] is not stored and so on (the pattern 102 repeats indefinitely).

Example 1 Let us consider the node sampling. The system generated by PRELUDE is S, V, R, C with:

1. S = {τ1 = (10, 2, 0, 10), τ2 = (30, 5, 0, 30), τ3 = (60, 30, 1, 60)} Let us for instance consider the words for vs. vs is stored in two different buffers b2,1 for τ2 → τ1 and b2,3 for τ2 → τ3. vs is produced once by τ2 while τ1 executes three times. vs is consumed with an fby , therefore there are two cells in the buffer and the initial value is in the first cell. During the period of τ2, τ2.0 stores the value in the second cell so that it does not interfere with τ1. Then, τ2.1 writes in cell 1 while τ1 reads three times the previous value in cell 2. This behaviour repeats indefinitely.

2. V = {i, o1, o2, v f , vs} such that τ1.in = {i, vs}, τ1.out = {o1, v f }, τ2.in = {v f }, τ2.out = {vs}, τ3.in = {vs, v f } and τ3.out = {o2}, 3. R = {(τ1, {(0, 0)}, τ2), (τ1, {(0, 0)}, τ3), (τ2, {(0, 0)}, τ3)} 4

SCHEDMCORE schedulability analyzer

The second step of the development consists in choosing a multiprocessor scheduling policy and verifying off-line the correctness of the schedule. Schedulability analysis with SCHEDMCORE relies on the lsmc_converter tool, which transforms the task set S, R produced by PRELUDE3 into a finite automaton describing all the possible executions. The automaton is then verified either with the UPPAAL model checker [START_REF] Behrmann | A tutorial on UPPAAL[END_REF] or with an ad hoc exhaustive search we developped in C. The lsmc_converter can also compute off-line optimal scheduling strategies (either static or dynamic priority based), which are correct by construction.

Handling precedence constraints

Let us first define formally the conditions that need to be verified. Let J denote the infinite set of jobs J = {τ i.k , 1 ≤ i ≤ n, k ∈ N}. Given a schedule, we define two functions s, e : J → N where s(τ i.k ) is the start time and e(τ i.k ) is the completion time of τ i.k in the considered schedule. We say that a dependent task set is schedulable under a given scheduling policy if the schedule produced by this policy respects all the constraints of the task set and all the job precedence constraints: Definition 3 Let S, R be a dependent task set. It is schedulable under a given scheduling policy if and only if:

∀τ i.k , e(τ i.k ) ≤ d i.k ∧ s(τ i.k ) ≥ o i.k ∀τ i.k → τ j.k ′ , e(τ i.k ) ≤ s(τ j.k ′ )
There are mainly two different approaches available for handling precedence constraints. The first approach relies on the use of (binary) semaphore synchronizations: a semaphore is allocated for each precedence and the destination task of the precedence constraint must wait for the source task of the precedence constraint to release the semaphore before it can start its execution. In the second approach, as shown in [START_REF] Chetto | Dynamic scheduling of real-time tasks under precedence constraints[END_REF][START_REF] Forget | Scheduling dependent periodic tasks without synchronization mechanisms[END_REF], if there is a precedence constraint from τ i to τ j , it is sufficient to ensure that τ j always starts after τ i and that τ i has a higher priority than τ j to respect the constraint. Unfortunately, this is difficult to apply to the multiprocessor case because tasks can execute simultaneously on different processors, regardless of their relative priorities. We would have to force τ j to be released after the worst response time of τ i , which is clearly sub optimal.

We therefore choose to ensure the precedence constraints by (1) using explicit synchronization mechanisms, there is simply no transition in our automaton model that enables a task to start executing before its predecessors (precedence constraints are taken as is) (2) forcing the execution of a task to take the wcet. Forcing the wcet prevents scheduling anomaly issues where execution times lower than the wcet lead to a system becoming unschedulable on-line while it has been proved schedulable off-line [START_REF] Graham | Bounds on multiprocessing timing anomalies[END_REF] (the off-line analysis is based on the wcet).

Modeling system execution

The scheduling of the system can be represented as a sequence of configurations.

At time t, the state of the system is represented by conf(t) = conf(τ 1 , t), . . . , conf(τ n , t) where for i ∈

[1, n], conf(τ i , t) = (T c i (t), C c i (t), O c i (t), D c i (t)). O c i (t) = max(O i -t, 0) is the time remaining until the first re- lease of τ i . If O c i (t)
> 0, all the other parameters are null. O c i (t) = 0 means that at time t the task has been released (at least once), in which case the other parameters have a value.

T c i (t) = T i + ((t -O i ) mod T i ) is the time remaining until the next release of τ i , D c i (t) = max(0, T c i (t) -(T i -D i ))
is the time remaining until the next deadline and C c i (t) is the remaining execution time for the current job of τ i . Since all the decisions are taken at integer dates, it is sufficient to describe a discrete model as shown in [START_REF] Guan | Schedulability analysis of global fixed-priority or edf multiprocessor scheduling with symbolic model-checking[END_REF]. The equivalent sequence of configurations is: m 1 represents the periodic extended precedence constraints for τ 1 . For instance, τ 1.0 can only start executing at date 1, which is the date of completion of its predecessor τ 0.0 . Note that at time 5, there are two overlapping configurations: one for the completion of τ 0.0 ((0, 0)) and one for the release of τ 0.1 ((5, 1)). This holds at any release.

Example 2 The real-time execution of the task set below on the left (O

i = 0 and T i = D i ),

with a fixed priority policy on two cores, can be represented as the Gantt diagram below on the right:

Analysis for existing sub-optimal online policies

There are mainly two approaches for scheduling a task set on a multiprocessor platform: a global strategy, where there is a single queue from which the tasks with the highest priority are selected to execute on any of the available processors, or a (semi) partitioning strategy where an algorithm allocates a distinct subset of tasks (or jobs) to a (subset of) processor(s) (leading to different queues for different processors and local scheduling).

The current implementation of SCHEDMCORE schedulability analyzer is targeted for global strategies but could easily be extended to (semi) partitioned policies. Currently, four policies are supported: (1) FP (fixed priority, tasks are assigned a specified priority off-line), ( 2) gEDF (multiprocessor extension of Earliest Deadline First, priority is given according to the next absolute deadline of each job), ( 3) gLLF (multiprocessor extension of Least Laxity First, priority is assigned regularly according to the remaining execution time), and ( 4) LLREF (Largest Local Remaining Execution First, mixing fair execution and local laxity).

The lsmc_converter takes a task set (in our case the model generated by PRELUDE), verifies several necessary schedulability conditions, such as checking that the load does not exceed the maximum load and models its execution with a given scheduling policy as a finite automaton [START_REF] Cordovilla | Multiprocessor schedulability analyser[END_REF] depicted below: ok ko update config() not sched?

The automaton works as follows:

• Variable list prio(t) contains the list of tasks ordered by decreasing priority at time t. Priorities are assigned according to the chosen policy, this is the only part that depends on the scheduling policy. The task τ i executes on a processor at time t if it is active, is not blocked by a precedence constraint and its position in list prio(t) is less than the number of processors;

• The precedence constraints are stored in a bidimensional array tabCurrent;

• The loop transition labelled update tab updates the configurations (T c τ (t), C c τ (t), O c τ (t), D c τ (t)), the array tabCurrent and the list list prio(t) after one unit of time elapses;

• The transition labelled not sched? can be fired only if the variable sched is equal to false, which only occurs when the task set is not schedulable (i.e. when some of the constraints of Definition 3 are violated).

The schedulability analysis then simply consists in proving that the transition not sched? is never fired. lsmc_converter generates a model to be verified either by the UPPAAL model checker or by an ad hoc exploration algorithm programmed in C. Since the execution is deterministic (when two tasks have the same priority, the scheduler always choose one task in a deterministic way to have the highest priority), it is sufficient to explore the execution until a sufficient point or a violation of the temporal constraints. The window is known to be [0, H] for synchronous constrained deadline task sets and for FP [START_REF] Cucu-Grosjean | Exact schedulability tests for real-time scheduling of periodic tasks on unrelated multiprocessor platforms[END_REF]. Otherwise, since the pattern of repetition of a schedule is a multiple of H, it is sufficient to explore until finding a repetitive configuration. The performance of this approach is relatively good compared to others, see §3.3 of [START_REF] Cordovilla | Multiprocessor schedulability analyser[END_REF] for detailed figures.

Optimal off-line schedule

A task set S, R is feasible for a given architecture if there exists a schedule that respects the temporal and the precedence constraints of every task (regardless of any specific policy). A scheduling algorithm is optimal [START_REF] Hong | On-line scheduling of realtime tasks[END_REF], with respect to an architecture and a class of policies (e.g. preemptive/non-preemptive, static/dynamic priority, etc.), if it can schedule all the feasible task sets. Being optimal on multiprocessor requires to be clairvoyant [START_REF] Hong | On-line scheduling of realtime tasks[END_REF] meaning that the scheduler must know a priori the future events, in particular, release times. Thus, on-line policies that only know of currently active tasks (which is the case of most of the on-line policies) cannot be optimal.

A widespread approach in critical embedded scheduling is to compute a feasible schedule off-line (also called sequencing). This provides several interesting properties:

(1) all the costs are assessable before run-time, for instance the number of preemptions, migrations and all context switches (2) no semaphores or synchronisation mechanisms are required since the dispatcher knows the predefined instants where tasks are resumed or suspended (3) off-line schedules are tractable and provide a very simple and efficient scheduler implementation (4) off-line scheduling is optimal.

The lsmc_converter is able to construct off-line schedules that can be run by the SCHEDMCORE runner. The construction also relies on the use of an automaton. This time, the objective is to find a trace reaching the state sched.

search sched ∀k ≤ min(m, nbtaskawake), ∀i1, . . . , ik, C c ij > 0

UPDATE STATE (i1, . . . , ik) t -max(Oi) = 0 mod H STATE SAVE schedulable ∧ t -max(Oi) = 0 mod H ∧ SAME PATTERN
The automaton works as follows:

• The transition UPDATE STATE updates the values of the tuples

(T c τ (t), C c τ (t), O c τ (t), D c τ (t))
after one unit of time. It is a non deterministic transition since we consider any combination of k active tasks, where k is the minimum between the number of processors m and the number of active tasks. These k tasks are assumed to access the processors at time t;

• The loop transition STATE SAVE may be fired nondeterministically at time t = max i≤n (O i ) mod (H) (remember that H is the hyperperiod). If it is fired, the current configuration is saved;

• The condition on the transition SAME PATTERN states that the task set must be schedulable and the backup is exactly equal to the current configuration. This means that we have found a repetitive schedule.

Optimal fixed priority assignment In [START_REF] Cucu | Feasibility intervals for multiprocessor fixed-priority scheduling of arbitrary deadline periodic systems[END_REF], Cucu and Goossens note that it is possible to determine an optimal fixed priority assignment by enumerating all the possible priority assignments and checking the feasibility of the resulting schedule on the schedulability interval defined in the paper. This can easily be encoded in the model, with minor modifications.

SCHEDMCORE execution environment

Once the program generated by PRELUDE has been proved schedulable, it can be executed on the multicore architecture target using the SCHEDMCORE execution environment.

SCHEDMCORE philosophy and objectives

The SCHEDMCORE runner, a.k.a. lsmc_run, requires several inputs:

• The specification of a task set: this can be done in two ways, either using a descriptive text file (SCHEDMCORE task file) or a dynamic library produced by the compilation of a PRELUDE-generated C file;

• The number of cores to be used;

• The scheduling policy: SCHEDMCORE currently offers four on-line policies (FP, GEDF, gLLF and LL-REF) and an off-line sequencer, which sequences the task set from a static trace. The user can also easily create his own policy and plug it into the platform.

A SCHEDMCORE scheduling policy is a dynamic library, which follows a specified API that allows it to be plugged into the runtime system. Adding a new policy simply consists in writting a C program that implements the API, compiling it as a dynamic library and specifying that the runner should be linked to this particular library. Thanks to the high level primitives provided by SCHEDMCORE, writting the policy itself is fairly simple. As an example, our implementations of the EDF and RM policies consist in a single C file, less than 50 code lines long.

We opted for a time-triggered platform, the pluggable scheduler wakes up at fixed time periods and forces tasks to consume all their WCET. This design makes the SCHEDMCORE runtime system simple, deterministic and impervious to scheduling anomalies [START_REF] Graham | Bounds on multiprocessing timing anomalies[END_REF].

The main target platform is Linux but it should be portable to any system that provides: a preemptive fixedpriority scheduler with at least 5 priority levels, some synchronization mechanisms (semaphores and/or barrier) and the ability to bind a thread to a particular processor (a.k.a. processor affinity). In this paper, we use a standard Linux system, with POSIX thread, SCHED_FIFO realtime scheduler and processor affinity support. The runtime is implemented as a user space library, which means that the Linux kernel does not have to be patched (as it can be with other solutions [START_REF] Faggioli | An edf scheduling class for the linux[END_REF]). Thus, the maintenance of SCHEDMCORE does not need to follow the evolution of the underlying operating system and it can be ported to another OS that offers the same primitives.

SCHEDMCORE architecture

The SCHEDMCORE runtime architecture is shown in figure 3. It consists of 3 management threads, which control the user real-time tasks that execute in separate threads. The scheduler skeleton is a POSIX thread, which manages the communication with the operating system and offers all necessary functions to manage task accesses to system resources. The execution of this entity is controlled by a semaphore. If the semaphore is unlocked the scheduler skeleton can preempt other tasks and performs all necessary actions for the task scheduling.

The schedule function of the pluggable scheduler implements one scheduling policy and takes scheduling decisions. It consists of a dynamic library that must provide 3 mandatory functions: initialize, schedule and finalize. The generic scheduler skeleton calls those functions at appropriate instant. Thus, systems instructions and scheduling instructions are separated which makes writing a new pluggable scheduler easier.

The time manager entity is woken up periodically at a fixed time interval by the operating system and it unlocks the scheduler skeleton semaphores to trigger the scheduling of the tasks. We chose to separate the time manager from the scheduler skeleton because we may decide (in further work) to release the scheduler at some other instants, or to use the usual event-based triggering.

The design of SCHEDMCORE requires 5 operating system preemptive priority levels, that define how each thread can preempt the others. The first (highest) is assigned to the time manager task, the second is assigned to generic scheduler task, the third is the priority corresponding to Execution task state, the fourth to Idle task state and the last one is the one used by the main program. The task states will be explained in the following section.

Execution with an on-line scheduler

The task model implemented in the SCHEDMCORE runner can be represented by the automaton of Figure 4.

A task is Active when it is not blocked by a precedence constraint and ready to run (released). This corresponds

(τi, t) = false if O c τ i (t) > 0 or C c τ i (t) = 0 or prec(τi, t) true otherwise where prec(τi, t) = true ⇔      ∃j, k, τ j.k → τ i.(t-O i )/T i ∧ (C c j > 0 ∧ k = (t -Oj)/Tj) ∨k > (t -Oj)/Tj
The condition prec derives from Definition 1 and the remark that at time t, task τ i executes the job τ i.(t-Oi)/Ti . An active task can either be in state execution, when it executes on one of the m CPU cores, or in idle state, when it is ready to run but all cores are already occupied.

Otherwise, the task is Inactive. An inactive task can either be in state waiting time, if the task has not started yet (not released) or has already completed its current job execution, or in state waiting dep, if the task is blocked by some precedence constraint.

The main work of a pluggable scheduler is to decide when to change the state of the concerned tasks, that it to say to define when the transitions of the task automata are fired. This is defined in the schedule function of the pluggable scheduler, which is called by the skeleton scheduler.

Therefore, schedule function computes the ordered queue which corresponds to list prio. We use an additional function pos prio which gives the position of an active task τ in list prio. The conditions associated with the different states are:

• Execution: active(τ, t) ∧ pos prio(τ, t) < m;

• Idle: active(τ, t) ∧ pos prio(τ, t) ≥ m;

• Waiting Time: not active(τ, t) ∧ not prec(τ, t);

• Waiting Dep: prec(τ, t)

Two operating system primitives are used to manage the task executions, semaphores and operating system priority. Semaphores are used to distinguish between active tasks (semaphore unlocked) and inactive tasks (semaphore locked). The m tasks with highest priority (according to scheduling policy) are assigned the highest (user task) priority of the OS and can go to state execution. The remaining tasks will go to or remain in idle state.

Execution with an off-line scheduler (dispatcher)

The schedule of the system can also be computed by an off-line scheduler (called a dispatcher). This mode of execution uses a trace table to decide at each instant which task must execute on which CPU core. The execution trace is not computed by the SCHEDMCORE runtime, it is instead computed off-line using the methodology described Section 4. The lsmc_tracer tool takes an UP-PAAL counter-example trace as input and generates a C file which contains the corresponding sequencing information. This file is compiled with a generic dispatcher C file to produce a dynamic library that conforms to the pluggable scheduler interface (as for on-line schedulers, a dispatcher is designed as a pluggable scheduler).

The automaton representing the execution of a task with such an off-line schedule is described Figure 5. A task can only be in two different states, execution or waiting. The execution table simply describes the transitions of the tasks at each instant. Transition (4) unlocks the task semaphore and the system priority of the task becomes that of executing tasks. Transition (2) locks the semaphore and the system priority of the task becomes that of waiting tasks. Transitions (1) and ( 3) basically have no effect. We illustrate our framework on the FAS presented in Introduction. It can be programmed with PRELUDE as described Figure 6. The compiler generates a task set made up of of 19 tasks and 26 precedence constraints (the inputs and outputs are translated into tasks).

WAITING EXECUTION (1) (2) (3) (4)
The results of the schedulability analysis of the task set, performed by SCHEDMCORE, are provided Figure 7. LLREF fails to schedule the task set because it tries to balance task execution without taking into account the urgency of a task that is due to its successors.

Conclusion

We described the design of an open end-to-end framework for the development of multi-periodic real-time systems targeting multicore architectures. We think that bridging the gaps between formal high-level view of a system and its real-time execution facilitates rapid prototyping and experimentation in the domain of multicore systems for time-critical applications. SCHEDMCORE aims at filling this gap using COTS operating system.

In the future, we will extend the PRELUDE language to integrate mode automata [START_REF] Colac ¸o | A conservative extension of synchronous data-flow with state machines[END_REF]. We will also provide automatic ways for aggregating nodes when there are too imported node Gyro Acq ( gyro , t c : i n t ) r e t u r n s ( o : i n t ) wcet 1 0 ; imported node GPS Acq ( gps , t c : i n t ) r e t u r n s ( o : i n t ) wcet 1 0 ; imported node FDIR ( gyr , gps , s t r , gnc : i n t ) r e t u r n s ( t o p d e , t o g n c , t o t m : i n t ) wcet 2 0 ; imported node PDE ( f d i r , gnc : i n t ) r e t u r n s ( p d e o r d e r : i n t ) wcet 1 0 ; imported node GNC US ( f d i r , gyr , gps , s t r : i n t ) r e t u r n s ( o : i n t ) wcet 2 1 0 ; imported node GNC DS ( u s : i n t ) r e t u r n s ( pde , s g s , pws : i n t ) wcet 3 0 0 ; imported node PWS( gnc : i n t ) r e t u r n s ( p w s o r d e r : i n t ) wcet 2 0 ; imported node SGS ( gnc : i n t ) r e t u r n s ( s g s o r d e r : i n t ) wcet 2 0 ; imported node S t r A c q ( s t r , t c : i n t ) r e t u r n s ( o : i n t ) wcet 2 0 0 ; imported node TM TC ( f r o m g r , f d i r : i n t ) r e t u r n s ( cmd : i n t ) wcet 1 0 0 0 ; s e n s o r g y r o wcet 0 ; s e n s o r g p s wcet 0 ; s e n s o r s t r wcet 0 ; s e n s o r t c wcet 0 ; a c t u a t o r pde wcet 0 ; a c t u a t o r s g s wcet 0 ; a c t u a t o r gnc wcet 0 ; a c t u a t o r pws wcet 0 ; a c t u a t o r tm wcet 0 ;

node FAS ( g y r o : r a t e ( 1 0 0 , 0 ) ; g p s : r a t e ( 1 0 0 0 , 0 ) ; s t r : r a t e ( 1 0 0 0 0 , 0 ) ; t c : r a t e ( 1 0 0 0 0 , 0 ) ) r e t u r n s ( pde , s g s ; gnc : due 3 0 0 ; pws , tm ) var g y r o a c q , g p s a c q , s t r a c q , f d i r p d e , f d i r g n c , f d i r t m , g n c p d e , g n c s g s , gnc pws ; l e t g y r o a c q = Gyro Acq ( gyro , ( 0 f b y tm ) * ˆ1 0 0 ) ; g p s a c q = GPS Acq ( gps , ( 0 f b y tm ) * ˆ1 0 ) ; s t r a c q = S t r A c q ( s t r , 0 f b y tm ) ; ( f d i r p d e , f d i r g n c , f d i r t m ) = FDIR ( g y r o a c q , g p s a c q * ˆ1 0 , ( 0 f b y s t r a c q ) * ˆ1 0 0 , ( 0 f b y gnc ) * ˆ1 0 ) ; gnc =GNC US ( f d i r g n c / ˆ1 0 , g y r o a c q / ˆ1 0 , g p s a c q , s t r a c q * ˆ1 0 ) ; many. We will improve the performance of the off-line schedule generator and we will enrich the scheduling policy libraries.

Figure 2 . Development process 1 . 2 . 1

 2121 Figure 2. Development process

τ 1 τ 2 τ 3 .Figure 3 .

 233 Figure 3. SCHEDMCORE runtime architecture

Figure 4 .

 4 Figure 4. On-line task

Figure 5 .

 5 Figure 5. Off-line task

  ( g n c p d e , g n c s g s , gnc pws ) =GNC DS ( gnc ) ; pde = PDE ( f d i r p d e , ( 0 f b y g n c p d e ) * ˆ1 0 ) ; s g s = SGS ( g n c s g s ) ; pws=PWS( gnc pws ˜> 1 / 2 ) ; tm = TM TC ( t c , f d i r t m / ˆ1 0 0 ) ; t e l

Figure

  Figure 6. The FAS in PRELUDE

Figure 7 .

 7 Figure 7. Schedulability analysis of the FAS

  is the (partial) communication function, where W is the set of dependence words.Periodic tasks A periodic task set S = {τ i } i=1,...,n is a finite set of tasks, where each task τ i has four real-time attributes (T i , C i , O i , D i ). T i is the period of repetition, C i is the worst case execution time estimation (wcet), O i is the first arrival time and D i is the relative deadline, with D i ≤ T i. We denote τ i.k the k th instance of τ

i (starting with instance 0), which we will call a job (or task job). The job τ i.k is released at date o i.k = O i + kT i . The hyper-period is the least common multiple of all the task periods, i.e. H = lcm 1≤i≤n (T i

  . C is defined as follows:

	i	o1	o2 vf	vs
	τ1 (1) (1)	b1,2 : (100)	b2,1 : (111222)
			b1,3 : (100000)
	τ2		b1,2 : (1)	b2,1 : (21)
				b2,3 : (10)

6. The FAS in PRELUDE

	policy	schedulable time in UPPAAL time in C
	FP (default) no	0.033s	0.000s
	gEDF	yes	45.734s	0.010s
	gLLF	yes	6.231s	0.029s
	LLREF	no	0.071s	0.000s
	OPT-FP	yes	out of memory	0.013s
	OPT	-	out of memory	out of memory

The PRELUDE compiler is available for download at http://www.lifl.fr/˜forget/prelude.html

The SCHEDMCORE environment is available for download at http://sites.onera.fr/schedmcore/

lsmc_converter accepts its own task file format description as well, one does not need PRELUDE in order to do schedulability analysis.


A Mezaouer
email: amezaouer@univ-blida.dz

J.-F Dupuy

K Boukhetala
email: kboukhetala@usthb.dz

A stratified log-rank test with missing stratum information and dependent censoring

The log-rank test is often used to compare randomized treatment groups with respect to the distribution of a failure time outcome. The so-called stratified log-rank test can be used when it is necessary to adjust for the effect of some discrete covariate that may be predictive of the outcome. In many applied situations, this discrete covariate is missing for some of the patients and moreover, the distribution of the censoring time depends on the treatment group. In this paper, we introduce a modified version of the stratified log-rank test, which accommodates both these problems simultaneously. The asymptotic distribution of this new test under the null hypothesis of equality of the randomized treatment groups is established. A numerical study is conducted to examine the finite-sample behavior of this test under both the null and alternative hypotheses.

Introduction

The log-rank test is widely used to compare randomized treatment groups with respect to the distribution of some failure time outcome. If one needs to control for a covariate that may be predictive of the outcome, and if this covariate is discrete, the log-rank can be generalized to the so-called stratified log-rank test (see, for example, [START_REF] Klein | Survival Analysis: Techniques for Censored and Truncated Data[END_REF] and [START_REF] Martinussen | Dynamic Regression Models for Survival Data[END_REF]). Consider a clinical trial where n patients are randomly assigned to K different treatment groups. We wish to compare survival between groups, while adjusting for some discrete factor S with L modalities (also called strata, such as income groups or disease stages for example). If λ k,l is the instantaneous hazard function for a patient in the kth treatment group and lth stratum, then the test for treatment effect can be formulated as H 0 : λ 1,l = . . . = λ K,l for every l = 1, . . . , L versus H a : "there exists j and j such that λ j,l = λ j ,l for some l". Let T 0 1 , . . . , T 0 n be the times from randomization to failure observed in the K pooled groups. Let C 1 , . . . , C n be right-censoring times (the C i are assumed to be independent of the T 0 i and non-informative). †Corresponding author

For each patient i, we observe T i = min(T 0 i , C i ) and ∆ i = 1(T 0 i ≤ C i ), where 1(•) is the indicator function. Assume that the data consist of n independent and identically distributed quadruplets (T i , ∆ i , G i , S i ), i = 1, . . . , n, where G i ∈ {1, . . . , K} and S i ∈ {1, . . . , L} respectively indicate the group and stratum of the ith patient. Let N i (t) = ∆ i 1(T i ≤ t), Y i (t) = 1(T i ≥ t), and define

E (n) k,l (t) = n i=1 Y i (t)1(G i = k)1(S i = l) n i=1 Y i (t)1(S i = l) .
Then the stratified log-rank statistic for the test of no randomized treatment effect is of the form U = (Z 1 , . . . , Z K-1 ) Θ-1 (Z 1 , . . . , Z K-1 ) ,

where for every k = 1, . . . , K -1,

Z k = n i=1 τ 0 1(G i = k) - L l=1 1(S i = l)E (n) k,l (t) dN i (t), (1)
τ denotes the end of the study period, and Θ is the estimated asymptotic covariance matrix of (Z 1 , . . . , Z K-1 ) . Under H 0 , U is asymptotically distributed as a χ 2 distribution with K -1 degrees of freedom (see, for example, [START_REF] Martinussen | Dynamic Regression Models for Survival Data[END_REF]).

In some applications, the stratum S may be missing for some patients. For example, consider the case where S represents the histological stage of patients included in a cancer clinical trial. The determination of S may require a biopsy, which due to expensiveness may not be performed on all the study subjects. One simple solution to handle such incomplete stratum information is to perform a complete-case analysis that is, to discard patients with unobserved stratum. This, however, may induce a substantial loss of power, as will be illustrated in our simulation study. Dupuy andLeconte (2008, 2009) considered the distinct but related problem of estimation in the stratified proportional hazards model with missing strata. The authors proposed a modified version of the maximum partial likelihood estimator, in which the unobserved stratum indicators are replaced by an estimate of their conditional expectation given available auxiliary covariates (this is the so-called regression calibration idea, see for example [START_REF] Carroll | Measurement Error in Nonlinear Models[END_REF]). Simulation results provided some evidence that such a replacement substantially improves on the complete-case analysis.

In many applications, it also happens that the distribution of censoring time depends on the treatment group. This arises, for example, when censoring follows from a study dropout caused by treatment toxicity. The treatment group with the heaviest toxicity will be more likely to have a higher dropout rate, and thus a higher censoring rate, than the other groups. Inverse probability of censoring weighted (IPCW) procedures have been proposed to remedy this problem (see for example [START_REF] Robins | Correcting for noncompliance and dependent censoring in an AIDS Clinical Trial with inverse probability of censoring weighted (IPCW) log-rank tests[END_REF], [START_REF] Yoshida | Estimation of treatment effect adjusting for dependent censoring using the IPCW method: an application to a large primary prevention study for coronary events (MEGA study)[END_REF], [START_REF] Cain | Inverse probability-of-censoring weights for the correction of time-varying noncompliance in the effect of randomized highly active antiretroviral therapy on incident AIDS or death[END_REF]).

In this paper, we propose and investigate a test of no randomized treatment effect, when the patients stratum information is only partially available and the distribution of censoring time depends on the treatment group. The test we propose combines the regression calibration and IPCW principles.

The rest of the article is organized as follows. In Section 2, we introduce the new test statistic and we derive its asymptotic distribution under H 0 . In Section 3, we describe a short simulation study investigating the finite-sample behaviour of the proposed test. A discussion concludes the paper in Section 4. An appendix contains the proofs of some intermediate technical results.

The test statistic

Assume that n independent patients are randomly assigned to K treatment groups. Assume that the stratum value is missing for some of these patients. Thus, a subsample is available where all variables (T, ∆, G, S) are observed, while only (T, ∆, G) are observed for the other patients. We assume that some auxiliary variables W ∈ R p are observed for all patients, and that W provides a partial information about S when S is missing. Let R be the indicator variable which is 1 if S is observed and 0 otherwise. Throughout the paper, we assume that T 0 and C are independent given G, S, W and R, and that C is independent of S and W given G. However as mentioned above, the distribution of the censoring time depends on the treatment group. We assume that T 0 is independent of W given S (that is, the auxiliary variables W provide no additional information about failure when the true stratum S is known), and that G is independent of S and W , as is the case in randomized clinical trials. We assume that R is independent of T 0 , C and G, and that 0 < P(R = 1) < 1. Finally, we assume that R and S are independent given W , which is the so-called missing-at-random assumption. In the sequel, this set of assumptions will be denoted by C1.

We consider the problem of implementing the stratified log-rank test of H 0 based on n independent vectors (T i , ∆ i , G i , W i , R i , R i S i ), i = 1, . . . , n of possibly incomplete data when moreover, the distribution of censoring time depends on the treatment group. To tackle simultaneously the missing strata and dependent censoring problems, we introduce a modified version of U , which is obtained by: (a) replacing any missing stratum indicator 1(S i = l) in (1) by its conditional expectation given the auxiliary W (this idea is related to regression calibration methods, see, for example, [START_REF] Thurston | Equivalence of regression calibration methods in main study/external validation study designs[END_REF], [START_REF] Weller | Regression calibration for logistic regression with multiple surrogates for one exposure[END_REF], [START_REF] Dupuy | A study of regression calibration in a partially observed stratified Cox model[END_REF]), and (b) weighting every patient by the inverse of the conditional (given the patient's treatment group) survival function of the censoring time (this idea is related to the inverse probability of censoring weighted principle, e.g., [START_REF] Robins | Correcting for noncompliance and dependent censoring in an AIDS Clinical Trial with inverse probability of censoring weighted (IPCW) log-rank tests[END_REF], [START_REF] Yoshida | Estimation of treatment effect adjusting for dependent censoring using the IPCW method: an application to a large primary prevention study for coronary events (MEGA study)[END_REF], [START_REF] Cain | Inverse probability-of-censoring weights for the correction of time-varying noncompliance in the effect of randomized highly active antiretroviral therapy on incident AIDS or death[END_REF]).

Precisely, we propose to base our test statistic on the following modified version of (1):

Z k = n i=1 τ 0 µ(G i , t) G k i - L l=1 D l i E (n) k,l (t) dN i (t), (2)
where for every i = 1, . . . , n, k = 1, . . . , K, l = 1, . . . , L, and t ∈ [0, τ ],

G k i = 1(G i = k), D l i = R i 1(S i = l) + (1 -R i )E[1(S i = l)|W i ], µ(G i , t) = 1/P(C ≥ t|G i ) (where P(C ≥ t|G) is the survival function of the censoring time in group G), and E (n) k,l (t) = S (n) k,l (t)/ S (n) l (t), with S (n) k,l (t) = 1 n n i=1 Y i (t)G k i D l i µ(G i , t) and S (n) l (t) = 1 n n i=1 Y i (t)D l i µ(G i , t).
Before stating our result, we need to introduce some further notations and some regularity conditions. For every i = 1, . . . , n and k = 1, . . . , K, let

V i,k = τ 0 µ(G i , t) G k i - L l=1 D l i E (n) k,l (t) dN i (t) - 1 n n j=1 L l=1 τ 0 Y i (t)D l i D l j µ(G i , t)µ(G j , t) S (n) l (t) G k i -E (n) k,l (t) dN j (t) (3)
and define

V i = (V i,1 , . . . , V i,K-1 ) and Σ = n i=1 V i V i . For every k = 1, . . . , K, l = 1, . . . , L, and t ∈ [0, τ ], define s k,l (t) = E[Y (t)G k D l µ(G, t)], s l (t) = E[Y (t)D l µ(G, t)],
and let e k,l (t) = s k,l (t)/ s l (t). The following regularity conditions will be needed in the proofs:

C2 There exists a positive constant c 0 such that P(C ≥ τ |G) > c 0 for every G ∈ {1, . . . , K}, and the survival function

t -→ P(C ≥ t|G) is continuous on [0, τ ]. C3 For every k = 1, . . . , K and l = 1, . . . , L, sup t∈[0,τ ] λ k,l (t) < c 1 for some finite positive constant c 1 .
C4 There exists a positive constant c 2 such that inf t∈[0,τ ] s l (t) > c 2 for every l = 1, . . . , L.

We are now in position to state the following result for the new test statistic

U := ( Z 1 , . . . , Z K-1 ) Σ -1 ( Z 1 , . . . , Z K-1 ) .
Theorem 1. Assume conditions C1-C4. Then under H 0 , as n → ∞, U converges in distribution to a χ 2 distribution with K -1 degrees of freedom.

Based on this theorem, the proposed test rejects

H 0 if U ≥ χ 2 1-α (K -1), where χ 2 1-α (K -1) is the quantile of order 1 -α of χ 2 (K -1).
Remark. The key to derive the null asymptotic distribution of the usual stratified log-rank test U is to represent Z k as a martingale process and to use a central limit theorem for martingales (see [START_REF] Fleming | Counting processes and survival analysis[END_REF] for example). This is not possible in our case, since the sum across strata in Z k has been replaced, in Z k , by a sum where each patient i such that R i = 0 contributes to each of the strata. Therefore, in order to obtain the null asymptotic distribution of U , we rather prove that Z k is asymptotically linear, and we use a central limit theorem for sums of i.i.d. terms.

Proof of Theorem 1.

Assume that all random variables are defined on a probability space (Ω, C, P), and define

F t,i = σ {N i (s), (1 -∆ i )1(T i ≤ s), G i , S i , W i : 0 ≤ s ≤ t}
as the σ-algebra generated by the event time and censoring histories of the ith patient over [0, t], and by the group, stratum, and auxiliary informations for this patient. Then by assumption, the F t,i -intensity of the counting process N i (t) is given by

Y i (t)λ i (t) = Y i (t) K k=1 L l=1 λ k,l (t)G k i 1(S i = l).
If S i is missing, the information for the ith patient is represented by the smaller σ-algebra

G t,i = σ {N i (s), (1 -∆ i )1(T i ≤ s), G i , W i : 0 ≤ s ≤ t} ⊆ F t,i . By the innovation theorem, the intensity of N i (t) with respect to G t,i is Y i (t)γ i (t) := E[Y i (t)λ i (t)|G t-,i ], where γ i (t) = K k=1 L l=1 λ k,l (t)G k i E[1(S i = l)|G t-,i ]. Finally, letting H t,i = (F t,i ) R i (G t,i ) 1-R i be the observed filtration, N i (t) has intensity Y i (t)ζ i (t) := Y i (t)[λ i (t)R i +γ i (t)(1-R i )] with respect to H t,i . It follows that M i (t) = N i (t) - t 0 Y i (s)ζ i (s)
ds is a martingale with respect to (H t,i ) t≥0 . In the sequel, we shall note κ l (t) = E[Y (t)ζ(t)D l µ(G, t)], l = 1, . . . , L (note that under the conditions stated above, κ l (t) < ∞ for every t).

The following lemma establishes a useful approximation of n -1 2 Z k . Its proof is given in the Appendix A. Lemma 1. For every i = 1, . . . , n and k = 1, . . . , K, let

Q i,k = τ 0 L l=1 D l i µ(G i , t) G k i -e k,l (t) dN i (t) -Y i (t) κ l (t) s l (t) dt . Then under C1-C4, n -1 2 Z k = n -1 2 n i=1 Q i,k + o p (1). Moreover, if H 0 holds, E[Q i,k ] = 0.
It follows from Lemma 1, from the multivariate central limit theorem, and Slutsky's theorem that under H 0 , n -1 2 ( Z 1 , . . . , Z K-1 ) converges in distribution, as n → ∞, to a (K -1)dimensional Gaussian vector with mean 0 and covariance matrix Σ

= E[Q 1 Q 1 ], where Q i = (Q i,1 , . . . , Q i,K-1 ) . Consequently, under H 0 , n -1 ( Z 1 , . . . , Z K-1 )Σ -1 ( Z 1 , . . . , Z K-1 ) d -→ χ 2 (K -1) as n → ∞. Σ however involves several unknown expectations. A consistent estimator of Σ is n -1 Σ := n -1 n i=1 V i V i , where V i = (V i,1 , . . . , V i,K-1 ) and V i,k is given by (3). To see this, since n -1 Σ = n -1 n i=1 (V i V i -Q i Q i )+n -1 n i=1 Q i Q i , it is sufficient to prove that V i,k -Q i,k p -→ 0 as n → ∞.
This follows from similar arguments and calculations as in the detailed proof of Lemma 2 (Appendix B). The details are therefore omitted.

Finally, it follows from Slutsky's theorem that

U := ( Z 1 , . . . , Z K-1 ) Σ -1 ( Z 1 , . . . , Z K-1 ) d -→ χ 2 (K -1) as n → ∞.

2

Remark. In practice, the weighting functions µ(G i , •) and/or the conditional probabilities E[1(S i = l)|W i ] may be either known from previous studies or completely unknown. In this latter case, estimated functions and probabilities have to be substituted in U . We come back to this issue in the simulation study.

A simulation study

We conducted a simulation study to evaluate the small to large-sample size behavior of the proposed test under various conditions. We considered the case of K = 2 randomized treatment groups and L = 2 strata. In each group and stratum, the event times T 0 i were generated from a Weibull distribution W (α, λ) with hazard rate λ(t) = αλt α-1 (the Weibull distribution is flexible and has a wide range of applications in survival analysis, see [START_REF] Klein | Survival Analysis: Techniques for Censored and Truncated Data[END_REF] for example). The failure times of stratum 1 in group 1 were generated from W (α 1 , λ 1 ), and those of stratum 1 in group 2 from W (α 1 , λ 1 r 1 ), where 'r 1 ' denotes the hazard rates ratio of two patients in stratum 1 of groups 1 and 2 respectively. The failure times of stratum 2 in group 1 were generated from W (α 2 , λ 2 ) and those of stratum 2 in group 2 from W (α 2 , λ 2 r 2 ), with 'r 2 ' being the hazard rates ratio of two patients in stratum 2 of groups 1 and 2 respectively.

We used α 1 = .5, α 2 = .75, λ 1 = .75, and λ 2 = 1.5. Three cases were considered for the pairs (r 1 , r 2 ) of hazard ratios: (a) (r 1 , r 2 ) = (1, 1), (b) (r 1 , r 2 ) = (1.5, 1.5), (c) (r 1 , r 2 ) = (1.25, 2). Case (a) corresponds to the null case of no difference between treatment groups, within each stratum. Cases (b) and (c) correspond to various magnitudes of difference between groups. In each case, the censoring times were generated from exponential distributions with parameters θ 1 in group 1 and θ 2 in group 2, with θ 1 and θ 2 chosen to yield censoring percentages equal to c 1 in group 1 and c 2 in group 2 (letting θ 1 = θ 2 ensures that the distribution of censoring depends on the treatment group).

Let n 1 and n 2 denote respectively the sample size in group 1 and 2 (with n = n 1 + n 2 ). We considered various values for (n 1 , n 2 ): (n 1 , n 2 ) = (50, 50), (n 1 , n 2 ) = (100, 100), and (n 1 , n 2 ) = (150, 150). The auxiliary variable W was taken to be 2-dimensional (W = (W 1 , W 2 ) ) with W 1 (respectively W 2 ) generated from the uniform distribution on [-1, 1] (respectively the normal distribution with mean 0 and standard deviation 0.5). A logistic regression model

P(S = 1|W ) = exp(b 0 + b 1 W 1 + b 2 W 2 2 ) 1 + exp(b 0 + b 1 W 1 + b 2 W 2 2 )
was taken for the relationship between S and W , with (b 0 , b 1 , b 2 ) chosen so that within each treatment group, each stratum contains approximately half of the patients.

The following stratum missingness percentages were considered: 20%, 40%. For each patient, the missingness indicator R was obtained by randomly drawing a Bernoulli random variable, with parameter chosen to yield the prescribed overall missingness percentage. The design parameters and their values are summarized in Table 1.

(S i = l|W i ) = E[1(S i = l)|W i ]
may either be known (from previous studies, for example), or completely unknown. In this latter case, one would estimate them and substitute the estimated values in the proposed test statistic U (the resulting statistic will be denoted by U in the sequel). However, the null asymptotic distribution of U may be somewhat distorted from the χ 2 (K -1) distribution. This, in turn, may affect the size and power of the test based on U . In fact, our simulation results show that as long as the µ(G i , •) and P(S i = l|W i ) are reasonably estimated, the prescribed level of the test is nearly maintained by U , and that U outperforms the complete-case log-rank test in term of power.

Various methods may be used to estimate P(S i = l|W i ). First, one may assume that there exists a known function h depending on an unknown parameter b, such that P(S i = l|W i ) = h(W i , b). Under the missing-at-random assumption, P(S = l|W, R = 0) = P(S = l|W, R = 1) = P(S = l|W ). One may thus obtain a consistent estimator b of b based on the patients {i : R i = 1} with observed stratum, and estimate h(W i , b) by h(W i , b). An alternative solution is, for example, to use the less restrictive local logistic regression approach (e.g., [START_REF] Loader | Local Regression and Likelihood[END_REF]).

Similarly, various approaches may be used to estimate µ(G i , t) = 1/P(C ≥ t|G i ). These include parametric and non-parametric methods.

The choice between competing estimation methods may be guided by the investigator's prior knowledge of the relationship between S and W , and of the censoring mechanism. In our simulation we considered the least favorable case where both P(S i = l|W i ) and µ(G i , •) are unknown. We estimated the P(S i = l|W i ) by local logistic regression, using the R package locfit (available at http://cran.r-project.org/web/packages/locfit/). We used non-parametric Kaplan-Meier estimators within each treatment group to estimate the censoring survival functions P(C ≥ t|G i ).

For each configuration of the design parameters, 1000 replications were obtained using the software R. Based on these 1000 repetitions, we obtained the empirical size (case (a)) and power (cases (b) and (c)) of the "estimated version" U of the proposed test U , at the significance level 0.05. For comparison, we included the results of the stratified log-rank test based on complete cases only (i.e. on individuals with known stratum). In the sequel, we shall refer this latter test to as U cc for short. Table 2 summarizes the results for an overall stratum missingness percentage equal to 40% (the results for 20% are similar and are therefore not presented). .5, 1.5) .823 .596 .762 .472 .822 .566 .698 .568 (1.25, 2) .772 .535 .813 .527 .746 .522 .611 .583 From these results, it appears that the proposed test U performs well and clearly outperforms the stratified test based on the complete cases (which was the only alternative test available so far for the problem considered in this paper). The empirical level of U tends to exceed (but only slightly) 0.05 in all cases. This may be due to the replacement of the unknown P(S i = l|W i ) and µ(G i , •) by their estimations, which causes the null asymptotic distribution of U to be slightly distorted from the χ 2 (K -1). But as expected, in cases (b) and (c), the powers of U are greater than those of U cc for every sample size and censoring percentages. In particular, U maintains a high power even when the censoring percentage heavily depends on the treatment group (20% in group 1, 50% in group 2), while at the same time the powers of U cc substantially decrease.

Discussion

We have constructed and investigated a modified version of the stratified log-rank test of no randomized treatment effect. This new test statistic is useful when the stratum information is missing at random for some patients and the distribution of the censoring time depends on the treatment group. From our simulations, we have found that this test performs well compared to the only alternative available so far, namely a complete-case based stratified log-rank test. Now, several questions still deserve attention.

First, we have assumed a missing-at-random mechanism for the stratum missingness. Investigating the robustness of the proposed test to a deviation to this assumption constitutes a topic for further numerical investigations. Extending the proposed method to non-ignorable missingness may be a non-trivial task however: the missing-at-random assumption is central in our proofs and for estimating the stratum belonging probabilities P(S i = l|W i ) when they are unknown. One may also investigate the case where the stratum missingness depends on the treatment group.

Second, in order to accommodate group-dependent censoring, we have used the inverse probability of censoring weighted principle, with weight function µ(G i , t) = h(P(C ≥ t|G i )) and h(x) = 1/x. Alternative test statistics Z k may be obtained by choosing other forms for h (we refer to [START_REF] Dirienzo | Bias Correction for Score Tests Arising from Misspecified Proportional Hazards Regression Models[END_REF] for alternative forms of weighting functions in the different context of bias correction for score tests arising from misspecified proportional hazards regression models). Searching for the function h which yields the most efficient testing procedure constitutes another non-trivial but very interesting task.

Another interesting future research direction is as follows. The stratified log-rank test can be viewed as a score test in a stratified Cox regression model. Therefore, the problem we considered in this paper can be viewed as the problem of implementing a score test of no randomized treatment effect in a stratified Cox model with missing stratum information. In the past two decades, a large amount of literature has been devoted to the problems of estimation and testing in the unstratified proportional hazards Cox model with missing covariates, but this literature has essentially focused on the case where the missing covariate has a proportional effect on the hazard of failure. As far as we know, the case where the missing covariate has a non-proportional effect and therefore is used to stratify was only few investigated. Dupuy andLeconte (2008, 2009) have considered this problem, but no treatment-dependent censoring was assumed. We hope that the present work can constitute a first step in the direction of relaxing this hypothesis.

Appendix A. Proof of Lemma 1.

We have

n -1 2 Z k = n -1 2 n i=1 τ 0 µ(G i , t)G k i dN i (t) -n -1 2 n i=1 L l=1 τ 0 µ(G i , t)D l i E (n) k,l (t) dN i (t) := A (n) 1,k -A (n) 2,k ,
and

A (n)
2,k can be written as

A (n) 2,k = n 1 2 L l=1 τ 0 e k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) -κ l (t)dt +n 1 2 L l=1 τ 0 v (n) k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) -κ l (t)dt + n 1 2 L l=1 τ 0 E (n) k,l (t)κ l (t)dt,
where v

(n) k,l (t) := E (n) k,l (t) -e k,l (t). We let the first term in A (n)
2,k unchanged. The second and third terms satisfy respectively the following two technical lemmas, whose proofs are postponed to the Appendix B:

Lemma 2. Under conditions C1-C4, n 1 2 L l=1 τ 0 v (n) k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) -κ l (t)dt p -→ 0 as n → ∞. Lemma 3. Under conditions C1-C4, as n → ∞, n 1 2 L l=1 τ 0 E (n) k,l (t)κ l (t)dt = n 1 2 L l=1 τ 0 e k,l (t) + S (n) k,l (t) s l (t) - s k,l (t) S (n) l (t) s l (t) 2 κ l (t)dt + o p (1).
Using these lemmas, we can re-write n -1 2 Z k as:

n -1 2 Z k = n -1 2 n i=1 τ 0 µ(G i , t)G k i dN i (t) -n 1 2 L l=1 τ 0 e k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) -κ l (t)dt -n 1 2 L l=1 τ 0 1 s l (t) s k,l (t) + S (n) k,l (t) - s k,l (t) S (n) l (t) s l (t) κ l (t)dt + o p (1).
Rearranging the terms in this expression concludes the proof of the first statement in Lemma 1. We now turn to the expectation of Q i,k under H 0 .

For every i = 1, . . . , n, k = 1, . . . , K, and

t ∈ [0, τ ], µ(G i , t)G k i is H t,i -measurable, hence the process ( t 0 µ(G i , s)G k i dM i (s)) t≥0 is a zero-mean martingale. It follows that E τ 0 L l=1 D l i µ(G i , t)G k i dN i (t) = E τ 0 µ(G i , t)G k i dN i (t) = τ 0 E µ(G i , t)G k i Y i (t)ζ i (t) dt,
where the first equality follows by noting that L l=1 D l i = 1. Similarly, we have

E τ 0 L l=1 D l i µ(G i , t) e k,l (t) dN i (t) = L l=1 τ 0 e k,l (t)κ l (t) dt, E τ 0 L l=1 D l i µ(G i , t)G k i Y i (t) κ l (t) s l (t) dt = L l=1 τ 0 e k,l (t)κ l (t) dt,
and

E τ 0 L l=1 D l i µ(G i , t) e k,l (t)Y i (t) κ l (t) s l (t) dt = L l=1 τ 0 e k,l (t)κ l (t) dt.
Thus

E[Q i,k ] = τ 0 E µ(G i , t)G k i Y i (t)ζ i (t) dt - L l=1 τ 0 e k,l (t)κ l (t) dt := B 1,k -B 2,k .
We now prove that under

H 0 , B 2,k = τ 0 E[G k ]E[Y (t)ζ(t)µ(G, t)] dt. First, remark that s k,l (t) = E Y (t)G k D l µ(G, t) = E E Y (t)G k D l µ(G, t)|G, S, W, R = E G k D l µ(G, t)E [Y (t)|G, S, W, R] H 0 = E G k D l µ(G, t)E 1(T 0 ≥ t)|S E [1(C ≥ t)|G]
where the third to last line follows from the assumptions C1 and from the fact that under H 0 , the distribution of T 0 does not depend on G. Then, by the independence of G and (S, W, R), we get that

s k,l (t) = E G k D l E 1(T 0 ≥ t)|S = E G k E D l E 1(T 0 ≥ t)|S .
Finally, using the same arguments as above, and the properties of conditional expectation, we have:

s k,l (t) = E G k E D l E 1(T 0 ≥ t)|S µ(G, t)E [1(C ≥ t)|G] = E G k E D l µ(G, t)E [Y (t)|G, S, W, R] = E G k s l (t). Thus, under H 0 , B 2,k = L l=1 τ 0 E G k κ l (t) dt = τ 0 E G k E[Y (t)ζ(t)µ(G, t)] dt.
Using similar arguments, we get that B

1,k = τ 0 E G k E[Y (t)ζ(t)µ(G, t)] dt under H 0 and thus, E[Q i,k ] = 0.

Appendix B. Proofs of Lemma 2 and Lemma 3.

Proof of Lemma 2. Let l ∈ {1, . . . , L}, and decompose

n 1 2 τ 0 v (n) k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) -κ l (t)dt as C (n) 1,k,l + C (n) 2,k,l , where C (n) 1,k,l := n 1 2 τ 0 v (n) k,l (t) 1 n n i=1 µ(G i , t)D l i dN i (t) - 1 n n i=1 µ(G i , t)D l i Y i (t)ζ i (t) dt and C (n) 2,k,l := τ 0 v (n) k,l (t) • n 1 2 1 n n i=1 µ(G i , t)D l i Y i (t)ζ i (t) dt -κ l (t)dt .
In the following development, we first show that C

(n) 1,k,l p -→ 0 as n → ∞. Note that C (n) 1,k,l is of the form: C (n) 1,k,l = n -1 2 n i=1 τ 0 H (n) i,k,l (t)dM i (t), where M i (t) = N i (t) - t 0 Y i (s)ζ i (s)ds and H (n) i,k,l (t) := v (n) k,l (t)µ(G i , t)D l i is a predictable process with respect to H t := n i=1 H t,i . Moreover, H (n) i,k,l (t) is bounded on [0, τ ] since |H (n) i,k,l (t)| ≤ |v (n) k,l (t)| • 1 c 0 ≤ 2 c 0 . Define the process (C (n) 1,k,l (t)) t≥0 by C (n) 1,k,l (t) = n -1 2 n i=1 t 0 H (n) i,k,l (s)dM i (s). Then (C (n) 1,k,l (t)) t≥0 is an H t -martingale, and
C (n) 1,k,l := C (n) 1,k,l (τ ). Now, the predictable variation process < C (n) 1,k,l > (t) of C (n) 1,k,l (t) is < C (n) 1,k,l > (t) = t 0 1 n n i=1 H (n) i,k,l (s) 2 Y i (s)ζ i (s)ds := t 0 X (n) (s)ds.
Under the regularity conditions stated in Section 2, it is not difficult to check that X

(n) (s) p -→ 0 as n → ∞ for every s ∈ [0, τ ]. Moreover, for any s ∈ [0, τ ] and n ≥ 1, |X (n) (s)| = 1 n n i=1 H (n) i,k,l (s) 2 Y i (s)ζ i (s) ≤ 1 n n i=1 H (n) i,k,l (s) 2 ζ i (s) ≤ 4 c 2 0 c 1
hence by Proposition II.5.3 of [START_REF] Andersen | Statistical models based on counting processes[END_REF], < C

(n) 1,k,l > (t) p -→ 0 as n → ∞. Next, for any > 0 and t ∈ [0, τ ], let C (n) 1,k,l, (t) = n -1 2 n i=1 t 0 H (n) i,k,l (s)1(n -1 2 |H (n) i,k,l (s)| ≥ )dM i (s).
Then the nonnegative < C

(n)

1,k,l, > (t) satisfies the following:

< C (n) 1,k,l, > (t) = t 0 1 n n i=1 H (n) i,k,l (s) 2 1(n -1 2 |H (n) i,k,l (s)| ≥ )Y i (s)ζ i (s)ds ≤ < C (n) 1,k,l > (t),
and therefore < C We now prove that C The second term on the right-hand side of ( 6) can be shown to converge to 0 in probability as n → ∞. The arguments are similar to the ones used in the proof of Lemma 2 and are therefore omitted.

  ,l, > (t) p -→ 0 as n → ∞. It follows (see Theorem 5.1.1 in Fleming and Harrington (1991)) that C (n) 1,k,l (t) d -→ 0 as n → ∞ for every t ∈ [0, τ ]. Therefore, C (n) 1,k,l (t) p -→ 0 as n → ∞ and in particular, if t = τ , C (n) 1,k,l p -→ 0 as n → ∞.

  n → ∞. To this end, we show that the integrand ofC i , t)D l i Y i (t)ζ i (t) -κ l (t) i , t)D l i Y i (t)ζ i (t) -κ l (t) . (4)Somewhat straightforward Glivenko-Cantelli arguments yield that E (n) k,l (t) is uniformly consistent (on [0, τ ]) for e k,l (t). It follows that as n → ∞, classes {Y (t), t ∈ [0, τ ]} and {µ(G, t), t ∈ [0, τ ]} are both Donsker (by Lemma 4.1 in Kosorok (2008)), and thus so is {µ(G, t)Y (t), t ∈ [0, τ ]} since products of bounded Donsker classes are Donsker. Similarly, {ζ(t), t ∈ [0, τ ]} is Donsker (Corollary 9.32 in Kosorok (2008)) and {D l } is Donsker. Finally, the class {µ(G, t)D l Y (t)ζ(t), t ∈ [0, τ ]} is Donsker (as a product of bounded Donsker classes) and therefore, the process G n(•) := n 1 2 n -1 n i=1 µ(G i , •)D l i Y i (•)ζ i (•) -κ l (•)converges weakly to some mean zero Gaussian process G. By the continuous mapping theorem, sup t∈[0,τ ] |G n (t)| converges weakly to sup t∈[0,τ ] |G(t)| and therefore sup t∈[0,τ ] |G n (t)| = O p (1). Combining this result and result (i , t)D l i Y i (t)ζ i (t) -κ l (t)

Table 1 .

 1 Design parameters and values included in the simulations

	Parameter Values	Description
	(r 1 , r 2 )	(1,1), (1.5,1.5), (1.25,2)	Ratios of hazard rates
	(c 1 , c 2 )	(5,20), (20,50), (30,20), (40,10) Censoring percentages
	n	100, 200, 300	Total sample size
	As mentioned previously, in practical situations the weighting functions µ(G i , •) and/or
	the conditional probabilities P

Table 2 .

 2 Empirical size and power of U and U cc , based on 1000 replicates.

					Censoring percentages (c 1 , c 2 )
			(5,20)	(20,50)	(30,20)	(40,10)
	n	(r 1 , r 2 )	U	U cc	U	U cc	U	U cc	U	U cc
	100 (1, 1)	.067 .045	.070 .067	.057 .045	.061 .057
		(1.5, 1.5) .354 .228	.272 .204	.401 .258	.360 .221
		(1.25, 2) .349 .227	.303 .220	.361 .231	.314 .216
	200 (1, 1)	.074 .055	.065 .051	.067 .049	.051 .047
		(1.5, 1.5) .692 .440	.565 .369	.651 .381	.574 .423
		(1.25, 2) .592 .386	.605 .418	.603 .403	.465 .401
	300 (1, 1)	.060 .046	.076 .048	.078 .052	.073 .048
		(1


Ferréol Soulez
email: ferreol.soulez@univ-lyon1.fr

Sébastien Bongard

Éric Thiébaut

Roland Bacon

RESTORATION OF HYPERSPECTRAL ASTRONOMICAL DATA FROM INTEGRAL FIELD SPECTROGRAPH

Keywords: Inverse problems, Deconvolution, Image restoration, ntegral field spectrograph, Astronomy

In this paper we present a method for hyper-spectral image restoration for integral field spectrographs (IFS) data. It takes advantage of all the spectral and spatial correlations in the observed scene to enhance the spatial resolution. We illustrate this method with simulations coming from the Multi Unit Spectroscopic Explorer (MUSE) instrument. It shows the clear increase of the spatial resolution provided by our method as well as its denoising capability.

INTRODUCTION

In the last decade, the integral field spectrographs (IFS) have become a popular tool for astronomical observation. Such instruments are now installed on all the main optical telescope facilities around the world. They provide spatially resolved spectra of a whole region of the sky, yielding (θ, λ) data cubes -with θ the 2D angular position and λ the wavelengththat have several hundreds of wavelength bins. With IFS, astronomical data enters the hyper-spectral era. Dedicated image reconstruction techniques are thus needed to take full advantage of the data gathered by these instruments. Because the light is spread out on multiple channels instead of being integrated on a single image, the information content is increased at the cost of a lower signal to noise or achievable resolution for the same exposure time. Furthermore, atmospheric turbulence and instrumental response often spatially blur the observations, degrading the spatial resolution.

First attempts to restore multi-channel images consisted in applying classical 2D restoration techniques like Wiener filter or Richardson-Lucy algorithm on each individual channel. The caveat of these approaches is to ignore the natural spectral correlations present in the data. The first restoration technique specifically dedicated to multichannel data [?] was a MMSE restoration filter based on the assumption that signal autocorrelation is spatially and spectrally separable. This assumption was later relaxed[?] and many other multichannel linear restoration filters have been proposed since. More recently, Fourier/Wavelet restoration techniques have been adapted to multispectral data[?, ?]. Combining demixing and restoration some authors[?, ?] achieve enhanced spatial resolution given the strong assumption that the observed scene is composed of only a few materials with unknown spectrum.

Most of the work on restoration of multi-spectral images is dedicated to remote sensing and color (RGB) images. Those methods can't easily be directly applied to astronomical data with its specific features like large dynamic range and strong sharp features (for example narrow emission lines or peaked sources). To the best of our knowledge, restoration techniques for multi-spectral astronomical images have only been proposed for (x, λ) data (slit spectrography) [?, ?] or (x, y, λ) data composed of slit spectrography scans [?]. At the current time, no similar techniques have been proposed for hyperspectral data observed simultaneously via an Integral Field Spectrograph (IFS).

We propose a new deconvolution method based on an inverse problem approach. It is very generic and exploits intrinsic continuities of hyper-spectral data. We suppose that a good estimation of the point spread function (PSF) is provided by other means (e.g. by calibration on the telescope guiding stars or on information from the adaptive optics system) and defer the blind deconvolution problem to a later time.

Our approach will be illustrated on data provided by the MUSE IFS simulator. Still in construction, the MUSE IFS[?] will be installed on the ESO Very Large Telescope (VLT). It is a "slicer" based IFS that covers in its wide field mode a 60 ′′ × 60 ′′ spectroscopic field-of-view subdivided into a grid of about 300 × 300 spatial elements (spaxels). To each spaxel corresponds a spectrum, obtained by dispersing the light on 3463 equally spaced spectral bins from 480 nm to 930 nm.

PROBLEM FORMULATION

We consider an observed data cube y of N λ monochromatic images of N Ω pixels. We model it from the parameters x following the equation:

y = H • x + e . (1)
with e the noise vector, and H the linear operator which approximates the convolution by the effective PSF and the sampling by the detector.

The PSF is assumed to be spatially shift invariant but its shape may vary with the wavelength. Similarly wavelengthwise PSF's may be centered at a location θ λ which depends on the wavelength so as to account for imperfect instrumental alignment and atmospheric differential refractive index (ADR). Furthermore PSF is not necessarily normalized in order to account for the variable throughput (atmospheric and instrumental transmission).

The parameters x describe the object 3-D distribution I obj (θ, λ) with a finite number of coefficients by means of expansion onto a basis of interpolation functions. The angular and spectral step sizes can be chosen to match the effective angular and spectral resolutions, in order to reduce the number of model parameters. We chose to control the effective number of free parameters by means of regularization and to take the same angular and spectral grid resolution than the data. The model parameters then simplify after discretization:

x k,ℓ = I obj (θ k , λ ℓ ) , (2)
where λ ℓ is the effective wavelength in the ℓ-th spectral channel and θ k is the k-th angular position in an evenly spaced rectangular grid of pixels.

MAXIMUM A POSTERIORI APPROACH

Deconvolution is a typical ill-posed problem [?] which can be solved by adding priors in a classical Maximum A Posteriori (MAP) approach. This is achieved by estimating the object x + that minimizes the cost function f (x):

x + = arg min x f (x) , (3)
f (x) = f data (x) + f prior (x) . (4)
This cost function f (x) is the sum of a likelihood penalty f data (x) ensuring the agreement between the model and the data y, and a regularization penalty f prior (x) introducing subjective a priori knowledge about the object.

Likelihood and Noise Statistics

Assuming Gaussian noise, the likelihood penalty reads:

f data (x) = [y -H • x] T • W err • [y -H • x] , (5)
where the weighting matrix W err = C -1 err is the inverse of the angular-spectral covariance of the noise and approximations.

Assuming uncorrelated noise, W err is diagonal and Eq. ( 5) simplifies to:

f data (x) = j,ℓ w j,ℓ [y -H • x] 2 j,ℓ
where 1/w j,ℓ is the noise variance of the measurements at pixel j and channel ℓ. This model can cope with nonstationary noise and can be used to express confidence on measurements on each pixel of the data. Since unmeasured data can be considered as having infinite variance, we can readily deal with missing or bad pixels as follows:

w j,ℓ def = Var(y j,ℓ ) -1 if y j,ℓ is measured, 0 otherwise. ( 6
)
This treatment of missing data is consistent because (i) it consistently accounts for unmeasured data and bad pixels, and (ii) it allows to properly expand the synthesized f.o.v. to avoid field aliasing and border artifacts caused by convolution using Fourier transform. Except for very low detector noise (< few e -per pixel), we can approximate the total noise (Gaussian detector noise plus Poisson noise) by a non stationary uncorrelated Gaussian noise [?]:

w j,ℓ def = γ max(y j,ℓ , 0) + σ 2 j,ℓ -1 if y j,ℓ is measured, 0 otherwise, (7)
where γ accounts for the quantization factor of the detector and σ 2 j,ℓ is the variance of other approximately Gaussian noise on the pixel (j, ℓ), like for example read-out noise.

Regularization

Astronomical data is mainly composed of bright objects (stars, galaxy) over a flat background. Most of the quite large MUSE field of view will thus contain only background. As a consequence, the data will be intrinsically spatially sparse. This spatial sparsity prior can be enforced by means of mixed norms[?, ?]:

f sparsity (x) = k   ℓ x 2 k,ℓ + ǫ 2 -ǫ   ( 8
)
where ǫ is a small real number (ǫ ≈ 10 -9 ) that ensures the derivability in 0 (hyperbolic approximation of the ℓ 1 norm). This regularization enforces spatial sparsity and spectral correlation since it favors when bright spaxels in each spectral channel are at the same spatial location. The regularization defined in Eq. ( 8) does not ensure the spectral continuity of the solution whereas in practice this solution should be relatively smooth along the spectral dimension excepted near emission and absorption lines. For that reason we introduce an additional regularization function:

f spectral (x) = k,ℓ (x k,ℓ -x k,ℓ-1 ) 2 + ζ 2 -ζ (9)
This regularization tends to smooth the spectra x k but preserve discontinuity where |x k,ℓx k,ℓ-1 | ≫ ζ. This situation is for example encountered at absorption or emission lines, which shall not be smoothed.

Owing to the large difference of dynamical range between spectral channel of astronomical images, these regularizations lead to over-regularization of bright features or under-regularization of fainter ones. For that reason, we rather suggest to apply these regularization functions to an spectrally whitened object x ′ :

x ′ k,ℓ = x k,ℓ /s ℓ ( 10
)
with s ℓ = x k,ℓ k the spatially averaged object spectrumk denotes averaging over pixel index k. To avoid dealing with non-linear regularization, we estimate the mean object spectrum directly from the data:

s ℓ = y j,ℓ j /η ℓ ( 11
)
with η ℓ = η(λ ℓ ) the effective throughput in ℓ-th spectral channel. Note that this approximation is justified because we do not attempt to perform spectral deconvolution and because we use the same wavelengths grid for the sought distribution and the data. Finally our regularization penalty writes:

f prior (x) = αf sparsity (x ′ ) + βf spectral (x ′ ) . ( 12
)
where α and β are the hyper-parameters that have to be tuned to set the importance of the priors.

ALGORITHM DESCRIPTION

In the presented experiments, we suppose that there is no cross-talk between spectral channels, the matrix H is block diagonal, and since the PSF is isoplanatic, applying H consists in N λ discrete spatial convolutions, one for each spectral channel. Due to the convolution process, flux from the object just outside of the field of view is measured inside data.

To take this fact correctly into account, the estimated object has to be spatially larger than the observed field of view. At least half of the PSF support must be added on each side of the observed field of view to form the restored field of view. Furthermore, as in practice the convolution is computed using FFT, the same field extrapolation prevents border artifacts due to the circulant approximation. In the end, the application of H requires N λ spatial FFTs. The hyper-parameters α and β are estimated by trial and error. The deconvolved data cube x + is the solution of Eq. ( 3). It requires the minimization of the cost function f (x) that involves a large number of parameters (> 1 × 10 6 ).

To that end, we use the VMLM-B algorithm[?] which is a limited memory variant of the variable metric method with BFGS updates[?]. This algorithm has proved its effectiveness for image reconstruction and only requires the computation of the penalty function being minimized together with its gradient. The memory requirement is a few times the size of the problem. The quality of the presented algorithm was assessed on data from the MUSE IFS simulator. This data is a part of 51 × 36 spaxels (pixels size: 0 ′ . ′ 2 × 0 ′ . ′ 2) of the whole MUSE data cube. It contains 3463 spectral channels comprised between 480 nm and 930 nm. The PSF, shown on Fig. 1(a) is a Gaussian with a full width at half maximum that varies from 0 ′ . ′ 75 (3.75 pixels) at the red end to 0 ′ . ′ 92 (4.6 pixels) at the blue end. In Fig. 1(b), the data summed all over the channels is displayed so that it can be compared with the ground truth in Fig. 1(d). The corresponding deconvolution presented in Fig. 1(c), clearly illustrates the gain of resolution provided by our method. Both the shapes of the central galaxy and of the one near the upper left corner are recovered. Figure ?? displays the spectra of the central spaxel of the galaxy from the data (blue), the deconvolution (dashed red) and the ground truth (black). Even though regularizations introduce some ex-pected bias, the restored spectra are close to the ground truth: most of the spectral features are preserved. We display in Figures ??(b-c) spectral cuts through the heart of the central galaxy materialized by the yellow line in Fig. 1(d). These figures show (θ, λ) images zoomed between 568 nm and 576 nm for the data, the restoration and the ground truth. These plots show the resolution gain provided by our algorithm: both brighter structures are separated, with the spectrum at θ = 43th column visible in the restoration that was not visible in the data. Futhermore, the noise observed around the galaxy has been drasticaly reduced by our method. (a) Spectra of the spaxel materialized by the green 'x' in Fig. 1(d) for the data (blue), the deconvolution (red) and the ground truth (black).

RESULTS

CONCLUSION

We present a new method for deconvolving hyperspectral data. By exploiting both spatial and spectral correlations present in the data, our method provides a strong spatial resolution enhancement and an effective denoising along the spectral dimension (given that we suppose the absence of channel cross-talk). Its deblurring performance is assessed on simulations. The results clearly demonstrate the advantages to process the whole data-cube instead of processing each wavelength channel independently. This increase of the spatial resolution would be of special interests for many astronomical subject such as the study of the kinematics of galaxies or weak lensing surveys, for which the true shape of the galaxy uncontaminated by the observational blur is of special interest. During this work, Ferréol Soulez has been funded by the French Agence Nationale de la Recherche on the MiTiV project (ANR-09-EMER-008-01).

Fig. 1 .

 1 Fig. 1. Images from the data, the deconvolution and the ground truth summed along the spectral dimension.

Fig. 2 .

 2 Fig. 2. (θ, λ) images of the cut materialized by the yellow line in Fig. 1(d) magnified between 568 nm and 576 nm.