Steven Varoumas

Tristan Crolard

WCET of OCaml Bytecode on Microcontrollers: An Automated Method and Its Formalisation

Keywords: 2012 ACM Subject Classification Computer systems organization → Embedded software; Computer systems organization → Real-time systems Worst-case execution time, microcontrollers, synchronous programming, bytecode, OCaml Digital Object Identifier 10.4230/OASIcs.WCET.2019.5

Considering the bytecode representation of a program written in a high-level programming language enables portability of its execution as well as a factorisation of various possible analyses of this program. In this article, we present a method for computing the worst-case execution time (WCET) of an embedded bytecode program fit to run on a microcontroller. Due to the simple memory model of such a device, this automated WCET computation relies only on a control-flow analysis of the program, and can be adapted to multiple models of microcontrollers. This method evaluates the bytecode program using concrete as well as partially unknown values, in order to estimate its longest execution time. We present a software tool, based on this method, that computes the WCET of a synchronous embedded OCaml program. One key contribution of this article is a mechanically checked formalisation of the aforementioned method over an idealised bytecode language, as well as its proof of correctness.

Introduction

Due to their low cost and efficient power use, microcontrollers are heavily used by the embedded system industry. Nonetheless, the very limited memory and power resources of these devices has lead developers to use traditional low-level languages such as C or assembly. While being precise and powerful, these languages often lack the hardware abstraction that would allow programmers that are not system programming experts to write various applications, from home automation to more critical applications. For this reason, many projects have emerged that make it possible to run, on microcontrollers with less than 10 kilo-bytes (kb) of RAM, programs written in higher-level languages such as Java [START_REF] Brouwers | Darjeeling, a Feature-rich VM for the Resource Poor[END_REF], Scheme [START_REF] St-Amour | PICOBIT: a compact Scheme system for microcontrollers[END_REF] [6], or OCaml [START_REF] Vaugon | Programming Microcontrollers in OCaml: The OCaPIC Project[END_REF]. The runtime environments developed in these projects usually include an optimised virtual machine that interprets bytecode of the language directly on the device. Beyond the usual hardware abstraction layer, these languages offer a higher-level expressiveness for the development of programs, and more guarantees over them (thanks to strong static type systems and automatic memory management among other features).

5:2 WCET of OCaml Bytecode on Microcontrollers

Timing constraints are frequent in embedded systems: they must often retain quick reaction times, in order to handle rapid changes in their environments. Computing the worst-case execution time (WCET) of the program's routine responsible for reacting to input stimuli is thus necessary, in order to ensure that no input is ignored during execution, and thus prevent a potential incorrect or dangerous behaviour from the system. A usual method of handling these timing constraints is to rely on schedulability analyses together with the use of real-time operating systems such as FreeRTOS [START_REF] Guan | Open Source FreeRTOS As a Case Study in Real-time Operating System Evolution[END_REF]. However, those systems are often quite heavy and not particularly adapted for microcontrollers with low memory resources. For this reason, the approach we present relies instead on the synchronous programming model. This model offers a lightweight reactive and concurrent programming model, particularly adapted to the applications of microcontrollers and their material limitations [START_REF] Varoumas | Concurrent Programming of Microcontrollers, a Virtual Machine Approach[END_REF], as well as advantageous properties for WCET computing: by construction, the generated code does not contain loops and does not need annotations.

In this article, we use a toolchain stemming from our previous works that consists in an OCaml virtual machine implementation fit for running on various models of microcontrollers, as well as a compatible synchronous programming extension to the OCaml language. This solution enjoys the high-level features of the OCaml programming language in order to develop safer programs for embedded system, and benefits from a model of concurrent programming adapted to the scarce memory resources of microcontrollers. Our approach, which uses a bytecode representation of embedded OCaml programs, provides portability, and enables a modularisation of various analyses done over these programs.

The main contribution presented in this article is an automated method that computes the WCET of a bytecode program on partially unknown inputs, together with its correctness proof. The method is implemented as a prototype tool called Bytecrawler which supports any OCaml bytecode instruction found in a compiled synchronous program. The formalisation is currently based on a small idealised bytecode language, all the proofs have however been mechanically checked.

Section 2 presents some preliminaries: we introduce the OMicroB virtual machine [START_REF] Varoumas | A Generic Virtual Machine Approach for Programming Microcontrollers: the OMicroB Project[END_REF] and OCaLustre, a synchronous extension to the OCaml language [START_REF] Varoumas | Concurrent Programming of Microcontrollers, a Virtual Machine Approach[END_REF]. Section 3 describes our method for computing the WCET of OCaml bytecode, and its implementation. Section 4 contains the formalisation and the correctness proof of the method. Section 5 concludes with a short discussion about related and future works.

Preliminaries

OMicroB: an OCaml virtual machine for microcontrollers

OMicroB [START_REF] Varoumas | A Generic Virtual Machine Approach for Programming Microcontrollers: the OMicroB Project[END_REF] is an implementation of the OCaml virtual machine (named ZAM -Zinc Abstract Machine [START_REF] Leroy | The ZINC experiment: an economical implementation of the ML language[END_REF]) dedicated to run OCaml bytecode on microcontrollers with very scarce resources (typically, less than 10kb of RAM and less than 100kb of flash memory).

Lightweight and configurable, it provides every feature of the language and its runtime, such as a strong static type system with type inference, the use of various programming paradigms (functional, imperative, modular, and object-oriented), as well as automatic memory management with garbage collection (GC). OMicroB is designed to be portable, due to its implementation in standard C, that we consider as a portable assembly language, since most of the target models of microcontrollers come with a dedicated C compiler. We were successful in porting OMicroB to AVR microcontrollers, and ports to ARM Cortex-M0 (used by BBC micro:bit cards) and PIC32 are promising works in progress.

The main component of OMicroB is an optimised stack-based bytecode interpreter, capable of running all the 149 bytecode instructions of the language. This interpreter handles various registers (such as a stack pointer, a code pointer, an accumulator, etc.) whose values are accessed and modified by the instructions of the running OCaml program. Moreover, the interpreter may invoke various low-level primitives (mostly dealing with input/output operations) that are directly implemented in C.

As shown in Figure 1, in order to execute an OCaml program on a microcontroller, the source code is first compiled into bytecode (by the standard ocamlc compiler), which is then "cleaned" by the ocamlclean tool which performs dead-code elimination. This bytecode program is then embedded into a C source code as an array of bytes by a tool called bc2c. It is then compiled together with the OCaml runtime library and interpreter by some platform-dependent C compiler and linker, and thus turned into an executable. For simulation purposes, the program can also be compiled for (and executed on) a generic PC. OMicroB offers a safer and more expressive model of programming than the classical C and assembly languages that are traditionally used to program microcontrollers. The robustness and the greater hardware abstractions gained by such a higher-level programming language makes it easier for developers that are not particularly well-versed in hardware specifics to develop programs for embedded systems.

OCaLustre: a synchronous extension of OCaml

Due to their constant interactions with their environment, embedded systems must react quickly to various stimuli. A model of concurrent programming, suitable for the size of the limited resources of microcontrollers might thus be a welcomed addition to the programming of microcontrollers [START_REF] Varoumas | Concurrent Programming of Microcontrollers, a Virtual Machine Approach[END_REF]. Therefore, we use OCaLustre [START_REF] Varoumas | OCaLustre : une extension synchrone d'OCaml pour la programmation de microcontrôleurs[END_REF], a synchronous programming extension of OCaml, inspired by the Lustre [START_REF] Caspi | Lustre: A Declarative Language for Programming Synchronous Systems[END_REF] language. Synchronous programming rely on the main principle that computations made during a logical instant are instantaneous: that is, the time required at each instant by the program to compute its output from its input should be ignored when reasoning about programs.

In OCaLustre, as in Lustre, the elementary synchronous component is called a node: a function that takes flows of values as input and computes output flows. The body of a node is a system of equations, defining the set of output or local values. For example, the count node of figure 2 takes an r flow as an input, and returns the cpt flow as its output. The -» operator, equivalent to the fby ("followed-by") Lustre operator, has the following semantics: x = a -» b means that x is equal to a for the first instant, and then to the previous value of b (i.e. the value computed for b at the preceding synchronous instant) for the subsequent instants. Therefore, count computes the series of natural numbers, reset to zero when the value of the r flow is true. During the compilation of an OCaLustre program, each node is separately turned into standard OCaml code, compatible with the ocamlc bytecode compiler, and thus with OMicroB. Our compilation method, derived from the single loop code generation used by most Lustre compilers, produces lightweight code, compatible with the resources of microcontrollers. It also provides guarantees over the program, such as the absence of causality loops (i.e. flows should not mutually depend on each other during the same instant).

W C E T

Worst case execution time analysis of OCaml bytecode

In order to compute the WCET of a synchronous OCaLustre program, we created Bytecrawler, a tool that computes an upper bound of the execution time of each synchronous instant. This tool has the main advantage of working over the bytecode instructions of the virtual machine, rather than the underlying native-code instructions. It relies on the fact that bytecode analyses and platform-dependent analyses are separated ( [START_REF] Harmon | A Survey of Worst-Case Execution Time Analysis for Real-Time Java[END_REF]) to free the application developer from needing to put additional annotations in their program: loop bounds for the program are not required due to the chosen programming model, and lower-level annotations, that can appear in the bytecode interpreter or the runtime library, are provided by the platform developer.

Bytecrawler WCET computation function

Values of variables that come from the environment of the program (such as values returned by electronic sensors) are unknown at compile-time, and thus introduce an external nondeterminism. Bytecrawler can thus be seen as an abstract bytecode interpreter, extended to deal with unknown values. It is based on the standard method of [START_REF] Puschner | Calculating the Maximum Execution Time of Real-Time Programs[END_REF] that considers all possible paths of the program to compute the costlier one, but Bytecrawler refines this evaluation by performing a hybrid execution that uses concrete values whenever possible. The program is normally executed by Bytecrawler with all values that are known at compiletime, but the evaluation function is generalised to handle unknown values: when reaching a conditional branching bytecode instruction (namely, BRANCHIF and BRANCHIFNOT) whose condition's value is unknown, Bytecrawler will explores the different possible paths. Figure 3 is an excerpt of Bytecrawler's WCET computation function. It computes an upper bound of the timing cost (in cycles) of a given OCaml program using a cost function that maps each bytecode instruction to its cycle count, and following the bytecode semantics. In particular, the CCALL instruction corresponds to the call to an I/O primitive (written in C): the return value of such a call is always unknown at compile time.

Note that, while this method might work for any OCaml program, some limitations quickly appear: any looping program whose number of iterations depends of an unknown value could not be bounded, and the triggering of the garbage collection (GC) algorithm could break the WCET estimation. For this reason, Bytecrawler is better suited to work only on the synchronous extension of OCaml, which enjoys the nice properties that instants do not contain any loop, and they handle only basic data types which never trigger the GC. Remark. Symbolic execution [START_REF] Ballabriga | Symbolic WCET computation[END_REF] might allow us to refine even more the WCET estimate by pruning unreachable paths. Such an analysis for a stack-based virtual machine seems more involved than for imperative languages. For instance, the symbolic execution of Java bytecode is clearly not straightforward [START_REF] Pasareanu | Symbolic PathFinder: symbolic execution of Java bytecode[END_REF]). Moreover, functional programs usually feature higher-order functions and the OCaml virtual machine is tailored for an efficient evaluation of such programs. A recent experiment with a symbolic execution for Haskell [START_REF] Hallahan | Building a Symbolic Execution Engine for Haskell[END_REF] shows that defunctionalisation into some intermediate first-order functional language might be required. While our synchronous extension does not currently include higher-order functions, translating back from this higher-order bytecode language to a first-order intermediate representation would already prove difficult. We thus keep the study of this method for a future work.

Cost function for bytecode instructions

The cost function which associates each bytecode instruction to a cycle count is a finite map, represented as a table. This table can be computed by a classic WCET tool (such as the Bound-T execution time analyser [START_REF] Holsti | Status of the Bound-T WCET tool[END_REF]), configured for the correct model of microcontroller.

A key benefit of our approach resides in the fact that this table must be computed only once for each microcontroller. The targeted devices are supposed to not induce timing anomalies [START_REF] Reineke | A Definition and Classification of Timing Anomalies[END_REF] due to their simple hardware model (cache-less, with in-order pipelines which are not timing anomalous [START_REF] Cassez | What is a Timing Anomaly?[END_REF]): this ensures that the computed cost of each bytecode instruction remains the same for every given OCaml program, and provides compositionality of the timing analysis.

For most of the bytecode instructions, the execution time is a constant, and our method thus does not overestimate their execution time. For some other instructions, the execution time is dependent on some argument (for example, the APPTERM instruction which performs tail calls depends on the number of parameters of the function) which is explicit in the bytecode: their cost can be pre-computed for every realistic instruction/parameter pair. For some of the remaining bytecode instructions (for instance, those dealing with higher-order closures), the execution time might be more difficult to bound statically. However, these instructions do not appear, by hypothesis, in the considered bytecode. Lastly, calls to C primitives (via the CCALL bytecode instruction) may require an execution time that depends W C E T 2 0 1 9 5:6

WCET of OCaml Bytecode on Microcontrollers

on the values of the arguments of the function. The WCET of these calls might thus overestimate the actual execution time of the CCALL evaluation. Handling these primitive calls relies on another cost table that maps each C primitive name to its maximum cycle count. The costs of these primitives are hardware-specific and need to be provided by the platform developer. Each primitive cost should be computed using classical methods and tools for WCET analysis.

Example. The excerpt of OCaml bytecode on the left side of figure 4 corresponds to the body of the step function created by compilation of the count node of figure 2. On the right side of this figure is displayed a table associating a cost to every bytecode instruction, which has been computed by Bound-T (configured for an AVR ATmega32U4 microcontroller). From this table, Bytecrawler computes a maximum cost of 2121 cycles for this function. For sake of simplicity, calls to I/O primitives (via the CCALL bytecode instruction) are not mentioned in this example. For instance, in a complete synchronous program, calls to input (resp. output) primitives happen in the beginning (resp. end) of each synchronous instant. These calls should thus also be taken into account.

Formalisation and correctness proof

In this section, we describe a formalisation of the inner working of the Bytecrawler tool. For this purpose, we first introduce an idealised bytecode language that serves as an illustration on a small subset of OCaml bytecode instructions. We then formalise how to compute an upper bound of the WCET for a program written with these instructions, and we prove the correctness of this computation. Our conjecture is that, because of the imperative nature of the OCaLustre programs, and because of the absence of dynamic memory allocation during the execution of a synchronous instant, the results over this idealised bytecode language can be transposed to the actual subset of the OCaml instructions coming from the compilation of an OCaLustre source code. The specifications presented in this section have been developed with the Ott tool [START_REF] Sewell | Effective tool support for the working semanticist[END_REF], and the various lemmas and theorems were proven using the Coq proof assistant [START_REF]The Coq Proof Assistant[END_REF]. The corresponding proof scripts are available at http://stevenvar.github.io.

The idealised bytecode language contains instructions for initialising variables (Init), assigning values to variables (Assign), arithmetic operations over integers (Sub and Add), branching (Branch and Branchif), as well as an instruction to terminate execution (Stop).

instr ::= Init x v | Assign x y | Add x y | Sub x y | Branch v | Branchif x v | Stop values, v, w ::= int_litteral | v + w | v -w
A state σ of a program is a tuple (P, pc, M ) that contains an array P representing the instructions of the program, a code pointer pc, and a memory M that associates any initialised variable name to its value. The small-step operational semantics of the language is defined on figure 5.

P[pc] = Init x v (P, pc, M ) -→ (P, pc + 1, M [x := v]) P[pc] = Assign x y M [y] = v (P, pc, M ) -→ (P, pc + 1, M [x := v]) P[pc] = Add x y M [x] = v M [y] = w (P, pc, M ) -→ (P, pc + 1, M [x := v+w]) P[pc] = Sub x y M [x] = v M [y] = w (P, pc, M ) -→ (P, pc + 1, M [x := v-w]) P[pc] = Branch v (P, pc, M ) -→ (P, v, M ) P[pc] = Branchif x v M [x] = 0 (P, pc, M ) -→ (P, pc + 1, M ) P[pc] = Branchif x v M [x] = 0 (P, pc, M ) -→ (P, v, M )
Figure 5 Operational semantics of the idealised bytecode language.

Definition 1 (Execution trace). An execution trace T is the sequence of states taken by the program, following the operational semantics of the language.

During the actual execution of a program, the values of all variables are known, and the corresponding execution trace is unique. Moreover, when the computation terminates normally, the trace ends on instruction Stop. We call such a trace deterministic: Definition 2 (Deterministic trace). The run function computes the deterministic execution trace beginning with a state σ:

P[pc] = Stop run((P, pc, M )) = [(P, pc, M )] σ -→ σ run(σ ) = T run(σ) = σ :: T Remark.
Note that, by definition, a deterministic trace is always finite. Indeed, since it is a certainty (due to absence of loops or recursion) that the actual execution of a synchronous instant will terminate, we are only interested in finite traces.

We associate a cost to every language instruction with the cost instr : instr → nat function. The cost of a transition is equal to the cost of the corresponding instruction, and the cost of a trace is the sum of the costs of all transitions:

cost step (P, pc, M ) = cost instr (P [pc]) cost(T ) = σ∈T cost step (σ)
W C E T 2 0 1 9

5:8

WCET of OCaml Bytecode on Microcontrollers

Erasure of variables

In order to compute an upper-bound for the execution time of an OCaLustre program, we reason over an abstract representation of the program which considers every possible path that the control flow can take. In this abstract representation, variables can hold unknown values, that represent the values that depend on the state of the program environment during execution. The grammar of the idealised bytecode language is extended accordingly, where unknown values are denoted by ⊥ :

values, v, w ::= int_litteral | v + w |v -w | ⊥
The memory of an abstract execution of a program might thus contain unknown values. We call such a memory an erasure of the actual memory of the program: Definition 3 (Erasure). A memory M is an erasure of a memory M (written M M ) if M and M share the same set of variables, but M contains more unknown variables than M . The following induction rules formally define the erasure of a memory:

∅ ∅ M M M ∪ x = v M ∪ x = v M M M ∪ x = v M ∪ x = ⊥
By extension, if two states σ and σ differ only by the fact that the memory of σ is an erasure of the memory of σ, we will use the same notation: σ σ and say that σ is an erasure of σ.

Since it depends only on the code pointer and the program instructions, the cost of a transition stays the same after an erasure of the memory of the program:

Lemma 4. ∀σ σ , σ σ ⇒ cost step (σ) = cost step (σ )
The wcet function of figure 3 computes the worst-case execution time of a program by always considering the maximum cost between two possible branches. In our formalism, this function is implemented by the cost max function as follows:

Definition 5 (Maximum cost). The cost max function computes the maximum cost of a program:

P[pc] = Stop cost max ((P, pc, M )) = cost instr (Stop) σ -→ σ cost step (σ) = c cost max (σ ) = k cost max (σ) = c + k P[pc] = Branchif x v M [x] = ⊥ cost step ((P, pc, M )) = c cost max ((P, pc + 1, M )) = k cost max ((P, v, M )) = k cost max ((P, pc, M )) = c + max(k, k )
The main result of this paper can now be stated as follows: the maximum cost of a program, computed from an erasure of its initial state, is greater than the actual execution time of the program.

Theorem 6 (Correctness). ∀ σ T , run(σ) = T ⇒ ∀ σ , σ σ ⇒ cost max (σ ) ≥ cost(T )
In particular, the cost max function can be applied to a memory where all variables are unknown (equivalent to the initial state of a synchronous instant) since this memory is an obvious erasure of all possible initial states. The result obtained is thus an upper bound of the execution times of all possible program runs.

The remainder of this section is devoted to the proof of this theorem.

Non-deterministic evaluation

The use of unknown values (⊥) induces a non-determinism on the evaluation semantics of the language: branching choices might depend on such variables, and thus it becomes statically impossible to guess which path would be the actual one. We thus introduce a non-deterministic operational semantics of the language which is a standard abstraction of the deterministic operational semantics where conditional branching is generalised to non-deterministic branching. The rules for this new semantics, denoted by , are the same as the rules for the deterministic semantics, extended with the two rules given in figure 6 (and the obvious lifting of arithmetic operations to account for ⊥).

P[pc] = Branchif x v M [x] = ⊥ (P, pc, M ) (P, pc + 1, M ) P[pc] = Branchif x v M [x] = ⊥ (P, pc, M ) (P, v, M )
Figure 6 Rules for non-deterministic branching.

Definition 7 (Non-deterministic trace). A non-deterministic trace is a finite sequence of transitions following the non-deterministic semantics of the language.

Because the non-deterministic semantics of the language is an extension of its deterministic semantics, a transition in the former stays the same in the latter. We say that it is embedded into the non-deterministic semantics:

Lemma 8. ∀ σ σ , (σ → σ ) ⇒ (σ σ )
Moreover, after erasing a memory, any transition from the deterministic semantics is still possible (possibly among other possible non-deterministic transitions):

Lemma 9. ∀σ 1 σ 2 σ 1 , (σ 1 σ 2 ) ∧ (σ 1 σ 1 ) ⇒ (∃ σ 2 , (σ 1 σ 2 ) ∧ (σ 2 σ 2 ))
By combining the two previous lemmas, we state that any erasure preserves the transition when going from the deterministic semantics to the non-deterministic one:

Theorem 10 (Preservation). ∀ σ 1 σ 2 σ 1 , (σ 1 → σ 2 ) ∧ (σ 1 σ 1 ) ⇒ (∃ σ 2 , (σ 1 σ 2 ) ∧ (σ 2 σ 2 ))
As illustrated by figure 7, the previous embedding and preservation properties hold for traces.

σ 1 σ 2 . . . σ n σ 1 σ 2 . . . σ n . . . σ 1 σ 2 σ n

Maximum trace

In order to prove that cost max is the cost of the "most expensive" run of the program, we need to formally define the notion of a maximum trace. We thus defined the run max function and proved that it actually computes a non-deterministic trace with the maximum cost.

Definition 11 (Maximum trace). The run max function is inductively defined as follows:

P[pc] = Stop run max ((P, pc, M )) = (P, pc, M ) :: ∅ P[pc] = Branchif x v M [x] = ⊥ run max ((P, pc + 1, M )) = T run max ((P, v, M )) = T cost(T ) ≤ cost(T ) run max ((P, pc, M )) = ((P, pc, M ) :: T ) σ -→ σ run max (σ ) = T run max (σ) = (σ :: T ) P[pc] = Branchif x v M [x] = ⊥ run max ((P, pc + 1, M )) = T run max ((P, v, M )) = T cost(T ) > cost(T ) run max ((P, pc, M )) = ((P, pc, M ) :: T )
As expected, cost max and run max are related by the following equality:

Lemma 12. ∀ σ, cost max (σ) = cost(run max (σ))
Using the trace preservation lemma, we derive the key property of this proof of correctness: the cost of the maximum trace is greater than the cost of every other finite trace, even when considering an erasure of their memories. In other words, the cost of the maximum trace computed from an erased state is greater than the actual cost of the program:

Lemma 13. ∀ σ T , run(σ) = (σ :: T ) ⇒ ∀ σ , σ σ ⇒ cost max (σ ) ≥ cost(σ :: T )
The combination of the previous lemmas allows us to conclude with our main theorem stating that the cost, estimated from an erasure of the initial memory, is an upper bound of the real cost of any actual run of the program: Theorem 14 (Correctness). ∀ σ T , run(σ) = T ⇒ ∀ σ , σ σ ⇒ cost max (σ ) ≥ cost(T )

Conclusion, related and future works

We presented a virtual machine approach to microcontrollers programming that makes it possible to run programs developed with a high-level, multi-paradigms, programming language on devices with limited resources. This language, OCaml, is extended with a synchronous programming model which is adapted to the concurrent nature of embedded programs. In particular, our approach provides a tool, called Bytecrawler, that makes it possible to compute the WCET of a synchronous instant by reasoning over the bytecode instructions of the program, thus providing an automated, factorised method to compute the reaction time of a program. The aforementioned method has been formalised on a idealised language, and a proof of its correctness has been developed in the Coq proof assistant. To the best of our knowledge, a formalised and automated WCET estimation of OCaml bytecode has not been done before. Similar studies exist on other embedded virtual machines for different programming languages (in particular, multiple works have been done on Java bytecode [START_REF] Bernat | Portable worst-case execution time analysis using Java Byte Code[END_REF][16] [START_REF] Harmon | A Survey of Worst-Case Execution Time Analysis for Real-Time Java[END_REF]), although these results do not seem to have been formally verified. Some recent work on proving the correctness of a WCET estimation tool is available [START_REF] Oliveira Maroneze | A Formally Verified WCET Estimation Tool[END_REF], but this tool is integrated within the CompCert C compiler, and thus operates only over C programs.

As a work in progress, we are adapting the formalisation and proof presented in this article to the actual, stack-based, language of OCaml bytecode instructions, in order to extract a certified implementation of Bytecrawler from Coq. As discussed in section 3, symbolic execution could improve the WCET estimate, but adapting it to the OCaml virtual machine might be difficult. Concolic execution [START_REF] Sen | Concolic testing[END_REF], a hybrid model of execution that mixes concrete and symbolic execution, seems close to our method, and might also refine WCET estimation.

Finally, it is worth mentioning that our solution could quite easily be adapted to other kinds of analysis: for example, modifying the cost function to indicate the number of words allocated in memory by each bytecode instruction could provide an upper bound for the total amount of memory used by the program, and thus be an estimate of the maximum memory usage of any given OCaml program.

Figure 1

 1 Figure1Compilation of an OCaml program with OMicroB (taken from[START_REF] Varoumas | A Generic Virtual Machine Approach for Programming Microcontrollers: the OMicroB Project[END_REF]).

Figure 2 A

 2 Figure 2 A synchronous node.

Figure 3

 3 Figure 3 Bytecrawler WCET computation function (excerpt).

Figure 7

 7 Figure 7 Trace preservation.


Sylvain Gerber
email: sylvain.gerber@mnhn.fr

USE AND MISUSE OF DISCRETE CHARACTER DATA FOR MORPHOSPACE AND DISPARITY ANALYSES

Keywords: character space, phylogenetic data matrices, disparity, macroevolution, missing data, morphological hypercube

The measurement of morphological variation in macroevolutionary studies is increasingly based on morphospaces constructed from discrete character data. This trend mostly results from the appropriation of phylogenetic data matrices as character spaces for carrying out disparity analyses. Phylogenetic matrices provide morphological descriptions of taxa as combinations of character states and thus appearif not conceptually, at least mathematically comparable to discrete character datasets found in numerical taxonomy or built for disparity purposes. Hence, phylogenetic matrices seem to constitute an abundant source of data readily available for morphospace analyses. Discrete character spaces have been generally described as more flexible than morphospaces capturing continuous shape variation. The discrete coding of morphology allows morphospaces to accommodate more disparate morphologies and the ability of discrete character schemes to handle missing data is also often emphasized. This flexibility comes at a cost, however. Multivariate ordinations of such spaces often provide deceptive visualizations and may invite the use of inappropriate methodologies for their exploration. The high amount of missing data that typifies many phylogenetic datasets is also problematic for the measurement of dissimilarity among taxa and can therefore be detrimental to the assessment of morphological disparity. In the present contribution, the properties of discrete character spaces are described and common pitfalls discussed. Graphical and methodological approaches are suggested to circumvent or limit their impact, and greater caution is recommended when using discrete character data for morphospace and disparity inferences.

The morphological disparity of a clade can be defined as a measure of the degree of morphological distinctness exhibited by its constitutive taxa (for reviews, see Foote 1997;[START_REF] Wills | Morphological disparity: a primer[END_REF][START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF][START_REF] Wagner | Evolutionary novelties[END_REF][START_REF] Hopkins | Morphological Disparity[END_REF]. Its study has been pioneered by palaeobiologists in the 90's (e.g. [START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF]Foote 1993), but the concept has since reached prominence as a macroevolutionary measure of morphological variation in both palaeontological and neontological contexts (e.g. [START_REF] Neige | Spatial patterns of disparity and diversity of the Recent cuttlefishes (Cephalopoda) across the Old World[END_REF][START_REF] Ricklefs | Cladogenesis and morphological diversification in passerine birds[END_REF][START_REF] Chartier | How (much) do flowers vary? Unbalanced disparity among flower functional modules and a mosaic pattern of morphospace occupation in the order Ericales[END_REF]. Methodologically, morphological disparity is estimated from the quantitative description of the distribution of taxa in morphospace, the space spanned by the set of morphological variables used to describe and compare the taxa studied.

Diverse kinds of morphospaces can be built depending on the descriptors used to characterise the sample of morphologies considered. In what follows, I focus on discrete character spaces.

These morphospaces have recently risen in popularity in palaeobiology, particularly among vertebrate specialists, as a consequence of the idea of co-opting phylogenetic data matrices as a source of morphological data for disparity analysis (e.g. [START_REF] Brusatte | Superiority, competition, and opportunism in the evolutionary radiation of dinosaurs[END_REF]Brusatte et al. , 2011;;[START_REF] Prentice | Evolution of morphological disparity in pterosaurs[END_REF][START_REF] Thorne | Resetting the evolution of marine reptiles at the Triassic-Jurassic boundary[END_REF][START_REF] Butler | How do geological sampling biases affect studies of morphological evolution in deep time? A case study of pterosaur (Reptilia: Archosauria) disparity[END_REF][START_REF] Ruta | The radiation of cynodonts and the ground plan of mammalian morphological diversity[END_REF]Hetterington et al. 2015;[START_REF] Romano | Disparity vs. diversity in Stegosauria (Dinosauria, Ornithischia): cranial and post-cranial sub-dataset provide different signals[END_REF]. [START_REF] Lloyd | Estimating morphological diversity and tempo with discrete charactertaxon matrices: implementation, challenges, progress, and future directions[END_REF] provided a recent account of the strengths and challenges associated with the use of discrete character data in studies of disparity and rates of morphological evolution. The theoretical and empirical features of discrete character spaces do not seem to have been much discussed however, and how these features affect morphospace and disparity studies has therefore been largely overlooked. In the present contribution, I describe and illustrate some properties of discrete character spaces and highlight potential pitfalls with the use of such spaces in disparity studies within current analytical frameworks. Although some of these issues do not seem to have ideal solutions for now, some suggestions are made to lower their impact on the interpretation of morphospace patterns and on estimates of morphological disparity.

The discussion mostly revolves around two points in need of further attention: the use of reduced-space ordination methods to visualize discrete character spaces and the large amount of missing data that typifies many phylogenetic datasets. Regarding the first point, it is often tempting to try to visualize the morphospaces we are working with. The use of ordination methods on discrete character spaces however leads to visualizations that can be prone to misinterpretations. It also tends to cause confusions between the morphospace itself and its ordination, inviting the use of inappropriate methods for analysing morphospace patterns. A hypercube graph example is developed to illustrate these aspects. The second point on missing data reflects a pressing concern given the widespread use of phylogenetic datasets for disparity analyses. Abundant missing entries have a substantial and detrimental effect on the assessment of morphological dissimilarity among specimens and can therefore strongly bias disparity estimates, in particular if those are calculated from ordinations of deficient pairwise distance matrices. Graphical options and weighting approaches are suggested to account for the uneven quality of specimen scoring.

DISCRETE CHARACTER DATA AND PHYLOGENETIC DATA MATRICES

Discrete character data have long been used as a means of quantifying morphological variation, including in the early years of the disparity research program (e.g. Foote 1992;[START_REF] Thomas | The skeleton space: a finite set of organic designs[END_REF]. Discrete coding corresponds to the morphological description of each specimen in a sample as a sequence of character states. The characters can possess two or more states, corresponding to binary and multistate characters respectively. In the latter case, the different states can be treated as ordered or unordered depending on the underlying mode of variation that is assumed.

One of the main advantages of discrete character schemes compared to other descriptors lies in their ability to accommodate a broader range of morphologies within a single analysis. Indeed, fairly disparate morphological structures, which can be challenging for morphometric approaches, can usually be more easily analysed within a discrete character framework.

Discrete coding also offers the possibility of including morphological features that are only present in a subset of the sample of specimens (a context that also precludes many alternative morphometric methods). The study of evolutionary novelties is an important research program in macroevolution (e.g. [START_REF] Wagner | Evolutionary novelties[END_REF]Erwin 2017) and the ability to account for the origination of novel morphological features (e.g. via a presence/absence coding scheme) is thus central to disparity. Discrete character data therefore appear especially relevant for studies at high taxonomic ranks where morphological variation is expected to be substantial and homology difficult to define uniformly across the taxa sampled. Discrete data are also useful if the morphological phenotypes considered include articulated or disconnected structures, such as the different bones of a skeleton. This flexibility generally comes at the expense of relatively cruder descriptions of morphologies, but careful coding schemes can still allow informative accounts of the morphological traits deemed relevant for a study.

Phylogenetic data matrices are instances of discrete character data and as such they have been increasingly considered as the basis for morphospace and disparity analyses in recent years.

One could nevertheless expect them to differ from the discrete character datasets that are specifically assembled for the purpose of disparity analyses. Phylogenetic data matrices are built for phylogenetic purposes and this has logical implications for the definition and coding of characters (e.g. [START_REF] Estabrook | An idealized concept of the true cladistic character[END_REF], notably regarding homoplasy (convergence and reversal) and autapomorphies (derived character state restricted to a unique terminal taxon).

The aim of phylogenetics is to establish branching orders and group taxa into evolutionarily coherent clades. Therefore, morphological features that would be known or anticipated to be highly homoplastic would not in principle be considered as suitable characters and included in a phylogenetic dataset. An opposite view on homoplasy is generally held in the field of disparity [START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF]Foote 1996). Far from being nuisances, convergences and reversals are seen as important evolutionary phenomena in studies of morphological disparity (see Foote 1996 for a discussion on total versus net changes) and their consideration, in conjunction with a phylogenetic hypothesis, improve the interpretations of evolutionary patterns. Autapomorphic characters are uninformative in parsimony analysis, which is still the most widespread approach in morphology-based phylogenetic inference (autapomorphies are nonetheless informative in maximum likelihood and Bayesian approaches; [START_REF] Lewis | A likelihood approach to estimating phylogeny from discrete morphological character data[END_REF]. One should therefore expect them to be rarely found in phylogenetic datasets. The view of their role in disparity studies ranges from central importance [START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF] to negligible [START_REF] Cisneros | Morphological diversity and biogeography of procolophonids (Amniota: Parareptilia)[END_REF]Ruta and Wills 2016) but for distinct conceptual and methodological reasons (see below).

Hence, phylogenetic and disparity datasets, respectively aiming at grouping and distinguishing taxa, in principle reflect two different modes of character matrix construction.

In a hypothetical phylogenetic dataset with no homoplasy and no autapomorphies, two terminal sister taxa (OTUs) share the same sequence of character states. They are, as scored, morphologically undistinguishable, which is not pertinent from a disparity viewpoint. The use of phylogenetic datasets for disparity purposes thus implicitly relies on the expectation of some erroneous assessments of character state homology and/or on the presence of parsimony-uninformative characters.

Empirically, it appears that discrete character datasets built for phylogenetic and disparity purposes do not in fact greatly differ from each other regarding the nature of the characters that they include. Based on measures of HER (Homoplasy Excess Ratio), [START_REF] Hughes | Clades reach highest morphological disparity early in their evolution[END_REF] showed that phylogenetic and disparity datasets on average exhibit comparable degrees of homoplasy, even though the small sample of datasets investigated did not allow further comparison of the statistical distributions of HER for the two kinds of datasets.

Studies of rates of character change in various groups have also shown that rate distributions are better fitted by gamma and lognormal models rather than by uniform (single-rate) models [START_REF] Wagner | Evolutionary patterns in early tetrapods. II. Differing constraints on available character space among clades[END_REF]Wagner 2012;[START_REF] Harrison | Among-character rate variation distributions in phylogenetic analysis of discrete morphological characters[END_REF]. These non-uniform distributions indicate significant rate heterogeneities among characters and suggest that phylogenetic data matrices commonly include homoplastic characters. Phylogenetic and disparity datasets also appear comparable with respect to the presence of autapomorphic characters. Most phylogenetic data matrices include autapomorphic characters indicating that phylogeneticists do not deliberately exclude them regardless of their relevance to the method of phylogenetic inference subsequently used (P. Wagner, pers. comm., based on a metaanalysis of 188 matrices of invertebrate clades showing that most matrices have at least one autapomorphic binary character, and nearly all have at least one autapomorphic character state on either a binary or multistate character).

Whether the two distinct research questions underlying the construction of phylogenetic and disparity datasets lead to comparable matrices or not, the use of the former as a surrogate for the latter still implies a series of reformatting steps. A phylogenetic data matrix cannot be directly co-opted without at least a few changes (see [START_REF] Hughes | Clades reach highest morphological disparity early in their evolution[END_REF]. For instance, the row corresponding to the outgroup needs to be removed in most cases. The question that motivated the construction of the phylogenetic data matrix might also have had an influence on the taxonomic coverage and (sub)sampling effort carried out, and it is important to assess the representativeness of the sampled OTUs for a disparity analysis meant to characterize the clade they belong to. Besides, it is also not uncommon for OTUs in a phylogenetic data matrix to represent different taxonomic ranks and some matrix rows might need to be amalgamated to level these differences and homogenise the dataset. Finally, phylogenetic data matrices conducive to phylogenetic trees with low level of homoplasy might indicate a stringent selection of characters and the possibility of expanding the character set in the prospect of a disparity analysis.

Hence, a phylogenetic data matrix can form a basis for the elaboration of a discrete character space, but the taxa and characters enclosed generally need critical re-evaluation. The characters eventually gathered generate a set of possibilities for morphological variation ("possibilities" in a combinatorial but not necessarily biological sense). This set can be finite or infinite depending on the definition of characters and their implied number of states (e.g., Hoyal Cuthill 2015). The next step is to define relationships among these possible morphological options, traditionally based on measures of morphological (dis)similarity, in order to generate a morphological space.

MEASURING DISSIMILARITY AND THE MORPHOLOGICAL HYPERCUBE

Hereafter, only binary characters are considered. This allows a simpler derivation of the observations and conclusions that are equally valid for geometrically more complex morphospaces built from datasets that include multistate characters. For convenience, it is also assumed that there is no missing data (this issue is dealt with in a following section) and that all characters can be treated equally. This last point means that all characters are given the same weight when estimating pairwise dissimilarity. Note that this equal weighting of characters is unlikely to be biologically meaningful given the anatomical diversity of characters typically combined, but this is the approach usually taken in empirical studies (see "CONCLUDING REMARKS AND PERSPECTIVES").

With discrete character data, the measurement of dissimilarity between two taxa is based on the number of matches and mismatches between their scored character states. For binary characters, this is equivalent to measuring the Hamming distance between the two sequences of character states. For instance, for a set of p=6 binary characters, the OTUs A (1, 0, 0, 1, 0, 1) andB (1, 1, 0, 0, 0, 1) have 4 matches (characters 1, 3, 5, and6) and2 mismatches (characters 2 and 4). The Hamming distance between A and B is thus 2. The distance between any possible pairs of morphologically distinct OTUs ranges from 1 to p. The Hamming distance satisfies the properties of equality, positivity, symmetry, and triangle inequality, and it is therefore a metric on the set of all sequences of binary character states of length p. The resulting discrete character space is a Hamming cube. If one considers the p states of all possible sequences as coordinates in p , one obtains a Euclidean representation of the discrete character space as a p-dimensional unit hypercube. A hypercube is the generalization of the concept of cube to higher dimensions. Each sequence occupies one of the vertices of the hypercube, and the distance between two sequences is the sum of unit-length edges separating their vertices (as opposed to the Euclidean distance between vertices). In this context, the morphospace is the spatial representation of taxa relative to an underlying set of possibilities for morphological variation. With p binary characters, this set of possibilities is the 2 p combinations of the character states corresponding to the 2 p vertices of the hypercube. Some geometric properties of such morphospaces can be easily inferred. For instance, adding an autapomorphic character to a dataset adds one dimension to the morphological hypercube, but it doubles the set of options for morphological variation: p+1 binary characters give 2 (p+1) sequences, that is, 22 p possible combinations of characters. The (p+1)-hypercube has now twice the number of vertices, but only one of these extra 2 p vertices is (and can be) occupied. This is not necessarily problematic but it must be kept in mind if one uses quantitative approaches whose performances are sensitive to morphospace dimensionality. Geometrically, an autapomorphy adds a dimension that separates the taxon that possesses it from all the other taxa of the clade. If autapomorphies are uniformly distributed across OTUs, their overall effect is to increase pairwise distances by pushing OTUs away from each other in orthogonal directions (in the Euclidean representation of the hypercube). In this case, their interest is limited for autapomorphies will not qualitatively affect the morphospace structure and relative measures of disparity extracted from it [START_REF] Cisneros | Morphological diversity and biogeography of procolophonids (Amniota: Parareptilia)[END_REF]. A more restricted and informative distribution of autapomorphies might also have a limited effect on overall disparity patterns inasmuch as their individual contributions to dissimilarity measures is relative to the total number of characters, but their interest rests upon their ability to single-out remarkable taxa and document the structure of morphospace occupation [START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF]Foote 1997).

Many other geometric features of discrete character spaces are not as easily grasped and highdimensional objects can be challenging to our intuition (e.g. Abbott 1884). Visualisation techniques such as projections onto lower-dimensional spaces are therefore frequently employed and it is important to assess the effect of such approaches on our understanding and interpretation of morphospace patterns.

MORPHOSPACE ORDINATION AND VISUALIZATION

Morphospace visualization is often considered as an important step for illustrating and discussing morphological patterns. Like most morphospaces however, discrete character spaces have a high number of dimensions and two-or three-dimensional representations of morphospaces consequently offer only a partial view of the amount and structure of morphospace occupation. Great care is therefore required when describing morphospace patterns from a reduced number of axes. This is especially true for discrete character spaces for which usually only a very small fraction of the variation is captured by the first axes of the multivariate ordination.

For discrete character data, it is customary to visualize a morphospace by carrying out a Principal Coordinate Analysis (PCoA) on the pairwise dissimilarity matrix D (an alternative not discussed here is Non-Metric Multidimensional Scaling, but the method aims at preserving ordering relations rather than distances among the objects it ordinates). PCoA provides a Euclidean representation of a dissimilarity matrix [START_REF] Gower | Some distance properties of latent root and vector methods used in multivariate analysis[END_REF], here analogous to the hypercube geometry described above. The principal coordinates generated (at most n-1

where n is the number of OTUs) are functions of the original variables (the morphological characters), as conveyed by the choice of the distance measure. If D is not Euclidean, because of the distance measure chosen (or because of missing data), negative eigenvalues will be produced. A transformation of D is therefore recommended prior to its ordination (e.g. square root transformation).

The resulting depiction of the morphospace is the set of taxa displayed as points in the space of the ordination, a generic representation common to all types of morphological descriptors.

Importantly however, this ordination is an embedding of the true discrete character space (morphological hypercube or similar but more complex discrete spaces) in a continuous Cartesian coordinate system, that of the principal coordinate axes. Unfortunately, this discrete structure is difficult to represent and is therefore often overlooked. Yet, conflating the geometric properties of the morphospace with those of its ordination has non-trivial consequences for the interpretation of morphological patterns and for the choice of methods employed to explore morphospaces.

A small example is presented in Figure 1 to illustrate these points. It is based on a subset of the archosaur dataset from [START_REF] Brusatte | Superiority, competition, and opportunism in the evolutionary radiation of dinosaurs[END_REF]. For the sake of the example, a stringent selection of characters had to be carried out in order to obtain a dataset where the OTUs retained are scored using only binary characters and without missing data. The culling of rows and columns from the original phylogenetic data matrix undoubtedly leads to an uninformative dataset, but it remains relevant for the present illustrative purpose. The subset considered includes 18 genera of crurotarsans coded for 10 binary characters (characters {2, 3, 4, 9, 19, 21, 22, 24, 25, 45} in Brusatte et al. 2008). As described above, dissimilarity among genera is calculated as the Hamming distance and these pairwise distances are stored in the distance matrix D. In the interest of visualization, a multivariate ordination of the morphospace is produced by carrying out a PCoA of the square root of D. Figure 1A shows the traditional representation of the morphospace with taxa abstracted as points in the principal coordinate space.

Two related questions are generally investigated or discussed when considering morphospace patterns. One concerns the overall amount of morphospace being occupied by a clade and the reasons for the existence of unoccupied regions. Discussions of this problem arose from Raup's seminal work on the unequal filling of the shell coiling morphospace [START_REF] Raup | Geometric analysis of shell coiling: general problems[END_REF] and are central to the field of theoretical morphology (e.g. [START_REF] Mcghee | The geometry of evolution: adaptive landscapes and theoretical morphospaces[END_REF]. The causes generally invoked include the functional and developmental limits imposed on the variational properties of phenotypes, ecological restrictions, and historical contingency (e.g. Raup 1972[START_REF] Maynard | Developmental constraints and evolution: a perspective from the Mountain Lake conference on development and evolution[END_REF][START_REF] Erwin | Disparity: morphological pattern and developmental context[END_REF][START_REF] Oyston | What limits the morphological disparity of clades?[END_REF]. The second research question aims at dissecting the temporal and/or spatial dynamics of morphospace occupation and forms the core of disparity analyses.

If the discrete character space is mistaken for its ordination, both these research questions can be misdirected. A priori unoccupied regions visible on the principal coordinate ordination might correspond to positions that do not belong to the discrete character space and therefore do not call for any biological explanation. Only a finite and non-randomly located (and in fact quite structured) positions correspond to the hypercube vertices, i.e. to meaningful positions in the actual morphospace. This is illustrated in Figure 1B, where the two-dimensional projection of the entire architecture of the 10-hypercube is displayed, that is, its 1024 vertices and its 5120 unit-length edges. To obtain this representation, all 1024 combinations of character states needs to be enumerated, along with a record of their 10 neighbouring combinations (those differing by only one character state change). With p characters, there are 2 p vertices and 2 p-1 p edges though, and it becomes rapidly difficult to enumerate all possible options and represent the hypercube graphically. Figure 1C shows a simplified representation where only a small subset of vertices and edges are computed and displayed, defined in a way that ensures a complete connectedness of the observed OTUs. It emphasizes the strong visual effect of the two-dimensional projection of a high-dimensional discrete morphospace. All the segments depicted represent edges and have therefore the same unit length despite their drastically variable projected lengths. Also, provided the distance matrix subjected to PCoA is Euclidean, edges are either parallel (if they correspond to the same character) or orthogonal, even if the ordination suggests otherwise. Spatial relationships among vertices inferred from the ordination should therefore be treated with caution since they will be misleading in most cases. This bias typically impacts graphical techniques such as convex hulls. The relative surface areas of these polygons can be poor indicators of the relative disparity of the groups they delineate in the ordination.

A possible way to visually circumscribe the morphospace within the space of principal coordinates could be to consider the hypersphere it is inscribed in (as conveyed by Fig. 1B).

By definition, all vertices are at a distance (p)/2 from the centroid of the hypercube. Hence, one could draw the circle of radius (p)/2 centred on the hypercube centroid to have an idea of the size of the morphospace. The centroid position can be calculated as the average principal coordinate scores of one of the observed vertices and of its diagonally opposite vertex (i.e. its opposite combination of character states). In practice however, the probability of having a good correlation between the true and the projected maximum distances is low (and decreases as p increases) and the circle drawn will generally largely exceeds the twodimensional distribution of the vertices. An alternative is proposed below based on the pattern of distribution of vertices.

Figure 1D illustrates the variation of the density of vertices under random orientations of the hypercube along the first principal coordinate axis (here the hypercube is centred at the origin of the principal coordinate system for convenience). It shows a density that decreases as one moves away from the hypercube centroid (this can also be noticed in Fig. 1B). This pattern has important implications: it means that even a uniform filling of the discrete morphospace will lead to a non-uniform distribution of taxa in the two-or three-dimensional visualization supplied by the PCoA (an orthogonal projection of the hypercube). Hence, density contour lines, heatmap representations, and similar approaches found in spatial point pattern analysis which might seem appropriate to document morphospace occupation (e.g. Ruta and Wills 2016) should be avoided or used with caution. More generally any method that relies on the continuous coordinate system of the ordination might be prone to similarly misleading description of morphological patterns.

The density of vertices along any principal coordinate axis is approximately distributed as N(0,0.25) (PCoA considers each binary character as a dimension of the Euclidean space, and their individual variance for the entire population of vertices is 0.25). Frequency ellipses can be used to delineate the bivariate distribution of the vertices in the plane of pairs of principal coordinate axes. In the bivariate case, a circular ellipse of radius ≈1.224 centred on the hypercube centroid will on average enclosed 95 percent of the morphospace vertices (with two degrees of freedom, (0.95) ≈ 5.991, so the radius of the circular ellipse is indeed (5.9910.25); see Fig. 1C). This sub-circular shape of the morphospace also explains the very low percentage of variance that is usually summarized by the first few principal coordinate axes. Hypercubes and related spaces do not show preferential directions of elongation (hence the smooth decrease of eigenvalues). Nevertheless, some datasets that are more structured and depart from a uniform filling of the hypercube can show a more efficient redistribution of variance following PCoA.

Such a structured morphospace occupation can be due to strong phylogenetic signal, clustered sets of character state combinations (morphologically distinct groups of taxa), or, more problematically, non-random distribution of missing data.

PHYLOMORPHOSPACE

phylogenetic datasets naturally encourage the use of their phylogenetic content in the course of disparity analyses, and the usefulness of contrasting patristic and phenetic dissimilarities in disparity studies has been emphasized (e.g. Foote 1996;[START_REF] Wagner | Patterns of morphologic diversification among the Rostroconchia[END_REF]). Here I restrict the discussion to the use of phylogenetic information to project the inferred tree onto the morphospace. This construction has been referred to as "phylomorphospace" [START_REF] Sidlauskas | Continuous and arrested morphological diversification in sister clades of characiform fishes: a phylomorphospace approach[END_REF], but see [START_REF] Strelin | The evolution of floral ontogenetic allometry in the Andean genus Caiophora (Loasaceae, subfam. Loasoideae)[END_REF] for a warning regarding the use of this term). It has been used to infer the placement of "ancestors" in the morphospace, assess the strength of the phylogenetic signal contained in a morphological dataset, and test for patterns of convergence (e.g. [START_REF] Klingenberg | Testing and quantifying phylogenetic signals and homoplasy in morphometric data[END_REF][START_REF] Stayton | The definition, recognition, and interpretation of convergent evolution, and two new measures for quantifying and assessing the significance of convergence[END_REF].

Producing a phylomorphospace visualization implies determining the location of the internal nodes of the tree in the morphospace, i.e., estimating their phenotypic trait values, and therefore requires a model of trait evolution. Importantly, the meaningfulness of these estimates depends on the adequacy of the model chosen. Different methods are employed for inferring ancestral trait values depending on whether the phenotypic traits considered show discrete or continuous modes of variation. In the case of discrete character spaces, it is therefore again critical to properly distinguish the discrete architecture of the morphospace from its ordination in a space equipped with a continuous coordinate system. Only methods devised for discrete traits (e.g. maximum parsimony or posterior probabilities calculated from the likelihoods of ancestral states in Mk models; [START_REF] Cunningham | Reconstructing ancestral character states: a critical reappraisal[END_REF][START_REF] Hunt | Integrating paleontological and phylogenetic approaches to macroevolution[END_REF] are appropriate for estimating ancestral character states in this context. Use of methods for continuous traits (e.g., squared-change parsimony, maximum likelihood models based on Brownian motion) applied to the principal coordinate scores of OTUs will lead to the incorrect placement of ancestral nodes in the ordination (e.g. [START_REF] Ruta | The radiation of cynodonts and the ground plan of mammalian morphological diversity[END_REF]Hopkins and Smith 2015). Their positions will be undefined locations in the discrete morphospace and subsequent analyses and interpretations relying on them will be flawed.

Figure 2 provides a simple theoretical illustration of this problem. Five OTUs (labelled A to E) are coded with 8 binary characters. These include 3 parsimony informative characters that support the phylogenetic reconstruction (Fig. 2A) and 5 autapomorphic characters that distinguish the OTUs from each other. Character states at internal nodes are estimated using maximum parsimony. As mentioned earlier, without autapomorphies and in the absence of homoplasy, sister OTUs would be merged with their supporting node. For instance, OTUs A and B and node i would all be described by the sequence (1, 0, 0).

Figure 2B shows the morphological hypercube (8 dimensions) embedded in the Euclidean space. The tips (OTUs) and the reconstructed internal nodes of the tree are located at their respective vertices of the hypercube. The branches of the tree correspond to edges (or sequences of edges in the general case) and are highlighted. This is a phylomorphospace representation that accounts for the discrete nature of the characters. Figure 2C shows the same hypercube and the same OTUs, but the distribution of internal nodes has been estimated using squared-change parsimony from the principal coordinate scores of the OTUs. The resulting depiction of the tree is substantially different from that of Figure 2B and none of the inferred internal nodes belong to vertices of the hypercube. This representation is invalid because it fails to acknowledge the true, underlying morphospace.

Phylomorphospaces built from continuous traits usually provide easily interpretable bivariate visualizations, with all internal nodes located within the convex hull defined by the OTUs and a degree of branch-crossing indicative of the amount of homoplasy. Conversely, with discrete characters, phylomorphospace visualizations tend to show very complex patterns from which little can be inferred (e.g. Halliday and Goswami 2015). Figure 2D illustrates the geometry of the phylomorphospace representation in the simple case considered. The parsimony informative characters supporting the tree define a subspace within the discrete morphospace.

It is also a hypercube but of lower dimensionality (here it is a cube). As described earlier and shown here, the autapomorphies push OTUs away from the subspace of parsimony informative characters in orthogonal directions. Brusatte et al. (2011) suggested the use of inferred ancestral nodes to derive phylogenetically corrected disparity curves. They described that the addition of reconstructed ancestors leads to an inflation of the morphospace, since these ancestors appear to fall outside the limits drawn by the observed taxa, as outlined by means of convex hulls. This is in fact a misinterpretation of the effect of the bivariate projection of the discrete morphospace and a drawback of the use of convex hulls in a discrete character context. These ancestors do not represent any aberrant or inexplicable combinations of character states but correspond to existing vertices on the discrete morphospace. Their apparent peripheral distribution results from the geometry of the morphospace (see Fig. 2B), the presence of homoplasy, and the distribution of the total set of taxa (observed and reconstructed) from which the ordination plane is defined.

THE PROBLEM OF MISSING DATA

Unscored entries in discrete character data matrices is problematic both for phylogenetic inference (e.g. [START_REF] Platnick | On missing entries in cladistic analysis[END_REF][START_REF] Maddison | Missing data versus missing characters in phylogenetic analysis[END_REF][START_REF] Wilkinson | Coping with abundant missing entries in phylogenetic inference using parsimony[END_REF]; but see [START_REF] Kearney | Fragmentary taxa, missing data, and ambiguity: mistaken assumptions and conclusions[END_REF] and for disparity analyses. They correspond to either missing or inapplicable data, and are traditionally reported as '?' and '-', respectively. Missing and inapplicable character states are handled similarly by phylogenetic software such as PAUP* [START_REF] Swofford | PAUP*. Phylogenetic analysis using parsimony (*and other methods)[END_REF] or TNT [START_REF] Goloboff | TNT, a free program for phylogenetic analysis[END_REF], and within standard disparity pipelines. The distinction between the two is not negligible from a morphospace perspective however.

Inapplicable character states

Within current methodological frameworks, it does not seem advisable to retain in a disparityaimed dataset characters that happen to be inapplicable for one or more taxa in a sample. These taxa do not have definable positions along the dimensions defined by the inapplicable characters but instead occupy the complementary subspace of the morphospace. A distance matrix can still be computed (albeit a deficient one) and ordinated by means of PCoA, but this will artificially enforce the placement of these taxa in a space they do not technically reside in. As it is sometimes done, disparity can even be calculated from this ordination, but the meaning of such measures is questionable on logical grounds. Novel kinds of morphospaces (and morphospace representations) need to be devised to accommodate the problem of inapplicable character states. This would be particularly useful for studies at high taxonomic ranks where OTUs tend to reside in only partially overlapping character spaces.

The effect of missing data

The central problem of missing data for disparity analysis is that it generates inconsistencies in the pairwise distance matrix D. If one considers the following three taxa with one missing entry for taxon A for the fifth characters:

Taxon A: (0, 0, 1, 0, ?) Taxon B: (0, 0, 1, 0, 0) Taxon C: (0, 0, 1, 0, 1) Comparison of sequences of B or C with A can only be done for four characters. From this partial information, C is seen as identical to A, B as identical to A, and yet, B and C differ by one mismatch. A taxon with missing data directly affects the n-1 distances it is involved in. It does not affect other pairwise comparisons, but it leads to inconsistencies in the overall distance matrix with cases as described above and a proportion of triplets of taxa whose pairwise comparisons do not satisfy the triangle inequality. If so, D is not fully metric anymore and its PCoA generates negative eigenvalues and potentially confusing visualizations (e.g. morphologically distinct taxa merged into a unique position). Some methods have been proposed to eliminate negative eigenvalues. These methods either add a specific constant value to the non-diagonal entries of D or D 2 (e.g. [START_REF] Lingoes | Some boundary conditions for a monotone analysis of symmetric matrices[END_REF][START_REF] Cailliez | The analytical solution of the additive constant problem[END_REF] or 'squeeze' the eigenvalues of the Gower-centred matrix (e.g., [START_REF] Hayes | Modification of estimates of parameters in the construction of genetic selection indices ('bending')[END_REF][START_REF] Kirkpatrick | Measuring selection and constraint in the evolution of growth[END_REF]. These corrections however differentially rescale dissimilarities depending on their magnitudes [START_REF] Legendre | Numerical Ecology[END_REF], and therefore need to be considered carefully in the context of disparity.

Obtaining a precise description of the effect of missing data on the distance matrix is a difficult task because numerous factors may affect, in interrelated ways, the response of the distance matrix to missing data. These factors include the shape of the frequency distribution of the pairwise distances (e.g., mean and variance), the shape of the character matrix M (n/p ratio), the distribution of missing data in M (amount and structure), and the use of scaling option for dissimilarities (when missing data occur, scaling is employed to make comparisons of dissimilarities fairer when these have been calculated from sequences of character states of different lengths. For instance, observed matches and mismatches between two taxa provide a dissimilarity measure which is then divided by the number of characters involved in its calculation -"mean character difference"; [START_REF] Sneath | Numerical Taxonomy[END_REF]Foote 1992).

The following simulation can help get a glimpse of the impact of missing data on the distance matrix. It is again a simplification of most real case studies, but it supplies an intuitive idea of the average behaviour of D while leaving aside the idiosyncrasies of empirical datasets. Let D be the Hamming distance matrix of a complete dataset (taken as a uniform distribution of n taxa in a p-hypercube), D k the distance matrix of the same dataset but plagued with a proportion k of uniformly distributed missing entries, and D Euc the Euclidean distance matrix calculated from the principal coordinates of D k . Results show that the correlations between D and D k , and between D and D Euc , both show a linear decrease with slope -1 as k increases (correlations measured as Pearson's r and Spearman's ). It means for instance that when k  10%, r and   0.9, and when k  50%, r and   0.5. Note that in this example, despite the changes in the ranking of distances in D k , morphological disparity measured as the mean pairwise distance (e.g. Foote 1992;[START_REF] Ciampaglio | Detecting changes in morphospace occupation patterns in the fossil record: characterization and analysis of measures of disparity[END_REF]) would on average be left unchanged. This does not occur for meaningful reasons however, but results from the uniform distribution of missing data. This random removal of entries in M has equal probabilities of increasing or decreasing a pairwise distance compared to its true value. As k increases, it also makes greater departures from the true value easier, leading to an increase in the variance of the distribution of distances (with the same mean) which explains the decrease in r and observed.

Patterns of missing data

Because of the non-random decay and fossilization of organismal characters (e.g. [START_REF] Sansom | Non-random decay of chordate characters causes bias in fossil interpretation[END_REF], missing data tend not to be randomly distributed in the character matrices but instead preferentially affect suites of characters that are generally lost together for anatomical, structural and/or compositional reasons. [START_REF] Smith | Joined at the hip: Linked characters and the problem of missing data in studies of disparity[END_REF] developed a linkage algorithm to account for the resulting differential quality of character scoring and showed that disparity estimates do indeed respond differently to random (Poisson distributed in their simulations) and linked missing data. In addition, taphonomic processes and preservational opportunities can also vary in time and space introducing additional heterogeneities in the pattern of missing entries among taxa.

All these factors can artificially inflate or reduce a measure of disparity depending on their particular combination and magnitude for the set of taxa considered. This potentially generates a collection of disparity measures (for instance corresponding to distinct time intervals, geographic areas or taxonomic groups) with reduced commensurability due to the limited overlap of their respective character subspaces. While a small amount of missing data does not seem to prevent reliable measures of disparity, it is unlikely to be the case when the proportion of missing data is significantly higher and heterogeneously distributed.

Unfortunately, the use of phylogenetic data matrices as disparity datasets has been accompanied by a substantial increase in the amount of missing data in disparity studies. [START_REF] Ciampaglio | Detecting changes in morphospace occupation patterns in the fossil record: characterization and analysis of measures of disparity[END_REF] provided a comparative survey of several disparity indices and investigated their sensitivity to various factors including missing data. At the time, they reported that disparity datasets were typically less than 25% incomplete and showed that the most widely employed disparity indices were relatively stable up to 30% of missing data.

With the incorporation of phylogenetic datasets within the pool of disparity datasets, this value of 25% is now regarded as an unrealistic upper bound [START_REF] Butler | How do geological sampling biases affect studies of morphological evolution in deep time? A case study of pterosaur (Reptilia: Archosauria) disparity[END_REF][START_REF] Lloyd | Estimating morphological diversity and tempo with discrete charactertaxon matrices: implementation, challenges, progress, and future directions[END_REF].

It is not rare to have phylogenetic datasets with more than 50% of missing data, with frequent instances of taxa with more than 80 or 90% of their entries unscored (e.g. [START_REF] Butler | How do geological sampling biases affect studies of morphological evolution in deep time? A case study of pterosaur (Reptilia: Archosauria) disparity[END_REF]).

Several approaches have therefore been recently suggested to cope with this intensified problem of missing data.

Suggested approaches

Phylogenetic correction methods to fill in missing data have been proposed and employed in empirical case studies (Brusatte et al. 2011;[START_REF] Butler | How do geological sampling biases affect studies of morphological evolution in deep time? A case study of pterosaur (Reptilia: Archosauria) disparity[END_REF]). One problem with such approaches is that they are dependent on the phylogenetic hypotheses and impose a model of character evolution on the reconstructed dataset. This might obscure the underlying dynamics of the true disparity signal and impede the test of macroevolutionary hypotheses. [START_REF] Smith | Joined at the hip: Linked characters and the problem of missing data in studies of disparity[END_REF] recommended that comparisons of disparity estimates among time intervals that involve sets of taxa differentially affected by missing data should account for this discrepancy. They suggested that taxa should be subjected to removal of some of their character entries, based on an adequate linkage pattern, so that all time intervals have as much missing data as the most incomplete one. While this rarefaction procedure might appear sensible, it is unlikely to improve the quality of disparity curves and may in fact decrease the signal-to-noise ratio. Missing entries are detrimental to the measurement of dissimilarity and increasing their amount is not a desirable option. If there are strong variations in the amount of missing data across time intervals, it is preferable not to artificially deteriorate reliable disparity estimates of well documented time intervals for the sake of comparing them to the low quality (and potentially meaningless) disparity estimates of other time intervals.

Recommendations

What other solutions then? Until new fossil discoveries help filling the gaps in the character matrix, there are options to minimize or at least provide reminders of the negative impact of missing data on the measurement of dissimilarities and the interpretation of morphospace and disparity patterns.

One approach to lower the relative amount of missing data is to work on subsets of characters.

If a pattern of character linkage applies to all the taxa investigated, then it means that the whole dataset can be partitioned into distinct subsets of characters with variable degrees of completeness. Reasonably complete subsets can then be combined or treated separately and their disparity patterns inferred. Beyond the preservational reasons underlying their definition, these subsets of characters might also reflect developmentally or functionally driven associations of traits into coherent units of evolutionary transformation (evolutionary modules; e.g. Schlosser and Wagner 2004). The disparity signals of these subsets might therefore be more helpful for identifying the determinants of a clade's disparity than its global disparity signal, because they are less likely to combine characters with different modes and rates of evolution [START_REF] Hopkins | Morphological Disparity[END_REF].

If one is interested in a global disparity signal, then approaches are needed to ensure that the missing data will not mislead the discussions of the results. Again, one risk comes from an undue reliance on morphospace ordination. Back to the hypercube example, if a taxon has k entries unscored in its p-length character state sequence, it means that the vertex the taxon truly occupies is one within a set of 2 k vertices, representing a proportion 2 k-p of the hypercube. As k increases, this can represent a very large set of possible locations (and a large proportion of the total set of vertices when k gets close to p). Not all possible options are equally likely of course, but still more than one can be expected to be plausible. An increasing proportion of such taxa with high k and high k/p lower the quality of the representation of D k in the principal coordinate space, introducing more and more pronounced discordances between disparity estimates measured from D k and from the ordination.

In general, despite the uneven quality and completeness of the description of taxa in the character matrix, an equal degree of confidence is implicitly attached to all taxa in the morphospace ordination: All of them contribute equally to the definition of the principal coordinate axes; all of them are depicted in the same way as unique points in the ordination space; all of them are considered equally informative in subsequent analyses carried out from the ordination scores. This uneven scoring quality needs to be documented and taken into account. A simple index q i =(p-k i )/p ranging from 0 to 1 (incomplete to complete) can provide a measure of the quality of the morphological description of taxon i.

For visualization purposes, it is possible to use a colour-code based on q (the column vector storing individual values q i ) in order to convey the degree of reliability one can attribute to a given taxon in the morphospace ordination (although its relative position can still be influenced by other taxa). It could also be possible to normalize q and use it as a vector of weights to calculate a weighted covariance matrix  of the principal coordinate scores. If there is substantial variation in scoring quality among taxa, a principal component analysis of  then allows a reorientation of the ordination space improving the representation of spatial relationships among the better-preserved taxa (with the reservations expressed earlier regarding the partial information provided by low-dimensional approximations).

For the construction of disparity curves, I do not recommend to carry out analyses from the morphospace ordination ("post-ordination disparity" in Lloyd 2016), but instead to obtain disparity measures directly from the distance matrix D k in order to minimize the impact of missing data. Each time interval is characterized by a set of taxa that is a subset of the total number of taxa sampled. The distance matrix of the subset is therefore a submatrix of D k .

Estimating disparity from this submatrix ensures that the estimate is unrelated to the amount and pattern of missing data of taxa occurring in other time intervals. Conversely, a disparity estimate inferred from the principal coordinate scores of the subset of taxa is influenced by all other taxa because the principal coordinate axes are defined from the whole D k . Hence, the detrimental effect of taxa with large amount of missing data will be spread out over the entire time range studied regardless of their temporal occurrence, and affect all ordination-based disparity estimates. Working with submatrices of D k restricts the impact of poorly-preserved taxa to their own time interval (the same argument holds if one is working with different taxonomic groups or geographic regions instead).

Disparity is conceived as the spacing and spread of taxa in morphospace. These features are traditionally calculated as the sum of variances and sum of ranges for a morphospace spanned by continuous variables. For discrete character spaces, the average pairwise distance is conceptually equivalent to the sum of variances (the latter can be shown to be proportional to the average squared pairwise Euclidean distance) and maximum pairwise distance indices analogous to range-based measures can be derived as well [START_REF] Hughes | Clades reach highest morphological disparity early in their evolution[END_REF]. The number of character state combinations is also an intuitive measure of the amount of morphospace occupied provided some of its weaknesses are accounted for (sensitivity to sample size and dimensionality. [START_REF] Thomas | The skeleton space: a finite set of organic designs[END_REF]Foote 1994;[START_REF] Ciampaglio | Detecting changes in morphospace occupation patterns in the fossil record: characterization and analysis of measures of disparity[END_REF].

Finally, it appears essential for disparity studies to systematically report the overall proportion of missing data and the average percentages of unscored entries per taxon and per character.

Maximum, minimum, and average values of q (or of any other equivalent quality index) should also be calculated for each time interval and plotted alongside disparity curves. This will provide a helpful indication to decide the level of confidence one can put in the disparity values reported.

CONCLUDING REMARKS AND PERSPECTIVES

On phylogenetic datasets as disparity datasets

Discrete character spaces have geometric features that make them quite distinct from other types of morphospaces encountered in evolutionary paleobiology. The discrete coding they allow is particularly relevant for the measurement of morphological variation at high taxonomic level. The co-option of phylogenetic data matrices as bases for morphospace construction has significantly increased the frequency of discrete character spaces in disparity studies. A positive outcome has been the broadening of the taxonomic scope of disparity analyses to include many vertebrate clades whose general morphology is difficult to characterize morphometrically. This rapid increase in the use of phylogenetic data matrices nevertheless calls for a critical assessment of their conceptual and operational value.

While it is often stated in these studies that "well established protocols" of disparity analyses were used, it is important to keep in mind that these protocols were devised for much larger datasets (in terms of number of OTUs). For example, the relative insensitivity to sample size of some disparity indices might need to be reconsidered if the number of taxa sampled in a time interval is very low (a frequent result of the temporal partition of small-size phylogenetic data matrices). Likewise, the often-mentioned ability of PCoA to handle missing data undoubtedly has limits, and obtaining principal coordinate scores does not guarantee their meaningfulness if the input distance matrix is poorly representative of the true pattern of variation. Hence, some aspects of the disparity pipeline need to be adjusted or reconsidered in order to account for the specificities of phylogenetic-based disparity datasets and ensure that their use is not associated with disparity inferences of lower quality.

On the measurement of dissimilarity

Despite the centrality of this problem (the choice of how to relate morphologies underlies the resulting topology of the morphospace), little has been said above on the possible options for measuring dissimilarity between pairs of OTUs. Only a few points are given below (and all assume the validity of a metric topology) but more work and methodological developments are needed in this area [START_REF] Gerber | Not all roads can be taken: development induces anisotropic accessibility in morphospace[END_REF](Gerber , 2017)). [START_REF] Lloyd | Estimating morphological diversity and tempo with discrete charactertaxon matrices: implementation, challenges, progress, and future directions[END_REF] proposed five criteria defining a good distance measure for disparity analyses:

(i) distances should be robust against missing data (this property can be assessed by simulations); (ii) the resulting distance matrix should be as Euclidean as possible (leading to no or few negative eigenvalues when carrying out PCoA); (iii) distances should be definable for all possible pairs of OTUs; (iv) their implied spatial pattern should be summarisable with as few dimensions as possible; (v) the set of distances should be as close to normally distributed as possible (because of the parametric nature of many statistical tests).

While points (i) and (iii) are legitimate, points (ii), (iv), and (v) are not justifiable. The choice of the distance measure should be biologically driven and not dictated by the requirements of ordination methods or of any subsequent analyses of the distance matrix. The notion of disparity emerges from the structure of the distribution of distances. Whether the latter departs from normality or does not lend itself to successful ordination in reduced-space should be a secondary concern (patterns of morphospace occupation can be inferred without visualisation). Defining a biologically meaningful measure of distance is a difficult task though given the anatomical diversity of the characters found in phylogenetic and disparity datasets. Weighting approaches might help account for the heterogeneity of characters and of their developmental bases (which invalidates equal weighting and impairs the value of simple measure such as the Hamming distance; see also [START_REF] Mitteroecker | The concept of morphospaces in evolutionary and developmental biology: mathematics and metaphors[END_REF]. The diverse developmental depth of characters was considered early on as an important issue calling for the differential weighting of characters when measuring dissimilarity between morphologies [START_REF] Gould | The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace[END_REF]. Equal weighting schemes seem to have been preferentially used so far, but treelike representations of development should in principle allow the definition of plausible hierarchies of characters and weighting schemes [START_REF] Schank | Generative Entrenchment and Evolution[END_REF][START_REF] Arthur | A theory of the evolution of development[END_REF]). The value of such hierarchical approaches needs to be evaluated.

On the visualization of morphospaces

It has been shown above why a morphospace must be distinguished from its ordination. In the case of discrete character spaces, the ordination consists in a Euclidean embedding of the morphospace and the ordination space is thus equipped with a coordinate system that differs from that of the morphospace. This distinction is important and both qualitative descriptions and quantitative analyses of the morphospace should not conflate the two spaces. Furthermore, in most empirical instances of discrete character spaces, the depiction of the spatial relationships among taxa supplied by the ordination is deceptive and its reliability can be further deteriorated by abundant missing data. Some graphical options have been suggested above to limit misinterpretations of morphospace patterns, such as the depiction of the connection structure of the vertices or the drawing of frequency ellipses of the distribution of vertices. It is also possible that the developmentally informed weighting schemes recommended above will improve the graphical representation of morphospaces by providing distributions of pairwise distances better summarized by standard ordination methods.

A morphological hypercube has been taken as an example throughout, but many datasets depart from this simple case. If multistate characters are present, the discrete space will be geometrically more complex than a hypercube and its architecture will depend on the ordered/unordered status of the character states. While this architecture might be more difficult to visualize, the options described above can be extended to these more complex cases (e.g. visualization of neighbouring vertices). Finally, obtaining informative visualizations of multidimensional spaces is a very general problem (e.g. [START_REF] Abramson | Two-dimensional projections of a hypercube[END_REF] and methods developed in different research fields might be fruitfully adapted to the context of morphospaces.

On missing data

Missing data appear to be the most serious problem brought on by the increasing popularity of discrete character data in disparity studies, even though its consequences are often overlooked. Numerous interacting biological and physical factors can be responsible for missing entries in discrete character matrices, making their specific effect on disparity estimates difficult to appreciate. Numerical or, ideally, analytical approaches are needed to better characterize the behaviour of dissimilarity and disparity measures in face of these perturbations. Progresses in this area could also help devise ways to lessen the impact of missing data. Meanwhile, reporting the structure (e.g. row and column distributions) and amount of missing data will allow a critical appraisal of the disparity patterns documented.

Combined with the small sample size issue resulting from the temporal partitioning of phylogenetic data matrices, missing data can lead to poorly supported disparity curves and temporal shifts in disparity between successive time intervals might not reflect any biological signal. In such contexts, coarser but more robust descriptions of clade shape seem preferable, for instance through the use of a centre-of-gravity approach (e.g. [START_REF] Foote | Morphological and taxonomic diversity in clade's history: the blastoid record and stochastic simulations[END_REF][START_REF] Hughes | Clades reach highest morphological disparity early in their evolution[END_REF]. While not as detailed as standard disparity curves, they will still contribute to our assessment of the relative frequencies of possible disparity patterns, and such assessments remain necessary for detecting potential trends and generalities in morphological macroevolution. B, the same ordination space but the true, discrete architecture of the morphospace, a 10hypercube, is shown by drawing the edges connecting neighbouring vertices. C, another visualization of the morphospace where only a subset of edges is represented. The subset is defined to ensure the complete connectedness of the occupied vertices. The size of points is proportional to the score on the third principal coordinate axis (orthogonal to the plane defined by the first two principal coordinate axes). The circular ellipse circumscribes ~95% of the locations of possible morphologies. D, mean variation of the density of vertices along a principal coordinate axis under 500 random orientations of the hypercube. The black line is the average of the 500 runs represented in grey.

FIGURE CAPTIONS

FIG. 2.

The mapping of phylogenetic tree onto morphospace. A, a simple phylogenetic hypothesis relating five taxa {A, B, C, D, E} and corresponding to the most parsimonious tree associated with the discrete character matrix shown on the right. The latter includes 3 parsimony informative characters and 5 autapomorphic characters. B, the visualization of the tree in the ordination space when the estimates of ancestral states at internal nodes {f, g, h, i} properly account for the discrete architecture of the morphospace, a 8-hypercube (inferred ancestral states are displayed in the matrix of the first panel). C, the visualization of the tree when the ancestral states are incorrectly inferred from the continuous Euclidean coordinate system of the ordination space (squared-change parsimony). D, same as B, but the subspace of parsimony informative characters is highlighted. It is a cube and autapomorphies push the five taxa away from it and from each other in orthogonal directions. Note that the projection of the morphospace here does not correspond to a principal coordinate analyses of the matrix of pairwise distances among {A, B, C, D, E}. It has instead been defined to better visualise the discrepancy between the two modes of ancestral states reconstruction, the architecture of the parsimony informative subspace and the effect of autapomorphic characters. Such an orientation could happen for instance with the presence of other groups of taxa and shows why reconstructed ancestors can appear "outside" the convex hull of the set of taxa from which ancestors are inferred.

FIG. 1 .

 1 FIG. 1. Morphospace and ordination in reduced-space. A, the traditional depiction of

  FIGURE 1 1055

Acknowledgements I am grateful to Jennifer Hoyal Cuthill and Gunther J. Eble for discussions on some of the topics covered in this contribution and to Zuzana Gerber and Yoland Savriama for comments on an earlier draft. Constructive comments by an anonymous reviewer and a detailed and insightful review by Peter J. Wagner helped improve and clarify the manuscript. Most of this work was done while holding a research associate position in the Department of Earth Sciences at the University of Cambridge and was funded by a Templeton World Charity Foundation grant (LBAG/143) awarded to Simon Conway Morris.