Bertrand Cloez

Benoîte De Saporta

Maud Joubaud

Optimal stopping for measure-valued piecewise deterministic Markov processes *

This paper investigates the random horizon optimal stopping problem for measurevalued piecewise deterministic Markov processes (PDMPs). This is motivated by population dynamics applications, when one wants to monitor some characteristics of the individuals in a small population. The population and its individual characteristics can be represented by a point measure. We first define a PDMP on a space of locally finite measures. Then we define a sequence of random horizon optimal stopping problems for such processes. We prove that the value function of the problems can be obtained by iterating some dynamic programming operator. Finally we prove on a simple counter-example that controlling the whole population is not equivalent to controlling a random lineage.

Introduction

Piecewise deterministic Markov processes (PDMPs) form a general class of non diffusion processes that was introduced by M. Davis in the 80's [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF][START_REF] Davis | Markov models and optimization[END_REF]. Such processes have deterministic trajectories punctuated by random jumps. They belong to the family of hybrid processes with a discrete component called mode or regime interacting with a Euclidean component. PDMPs can model a wide area of phenomena from insurance and queuing problems [START_REF] Davis | Markov models and optimization[END_REF], finance [START_REF] Bäuerle | Markov decision processes with applications to finance[END_REF], reliability [START_REF] De Saporta | Numerical methods for simulation and optimization of piecewise deterministic Markov processes: application to reliability[END_REF] to neuroscience [START_REF] Genadot | Multiscale piecewise deterministic Markov process in infinite dimension: central limit theorem and Langevin approximation[END_REF][START_REF] Riedler | Limit theorems for infinitedimensional piecewise deterministic Markov processes[END_REF], population dynamics [START_REF] Benaïm | Lotka-volterra with randomly fluctuating environments or "how switching between beneficial environments can make survival harder[END_REF][START_REF] Cloez | Probabilistic and piecewise deterministic models in biology[END_REF] and many other fields. In this paper we are especially interested in population dynamics applications. In this area, special cases of PDMPs include for instance growth-fragmentation processes for one or several interacting species [START_REF] Campillo | Links between deterministic and stochastic approaches for invasion in growth-fragmentation-death models[END_REF][START_REF] Cloez | Limit theorems for some branching measure-valued processes[END_REF][START_REF] Costa | A piecewise deterministic model for a prey-predator community[END_REF][START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF]. In that case, commonly, the deterministic part is the growth process that may depend e.g. on the age of the individual, on its size, on the quantity of available nutriment and the jumps correspond to fragmentation or division (for cells), birth or death events, abrupt changes in the environment,. . . In Davis' original construction, PDMPs are defined on subsets of R d , for some dimension d that may change when the process jumps. In this paper we are interested in extending the definition of PDMPs to measure-valued state spaces. Infinite dimensional PDMPs have already been introduced in [START_REF] Buckwar | An exact stochastic hybrid model of excitable membranes including spatio-temporal evolution[END_REF] (see also [START_REF] Genadot | Multiscale piecewise deterministic Markov process in infinite dimension: central limit theorem and Langevin approximation[END_REF][START_REF] Riedler | Limit theorems for infinitedimensional piecewise deterministic Markov processes[END_REF]). In those papers, PDMPs take values in a separable Hilbert space and model spatio-temporal phenomena occurring on neuronal membranes. Our approach differs as we are interested in measure-valued PDMPs to deal with population dynamics models. Instead of modeling the dynamics of only a single individual by a finite-dimensional PDMP, we aim at taking into account simultaneously the dynamics of all the individuals in the branching population when the population remains small and the stochastic approach is relevant and large scale approximations do not hold. Such a population can be represented by a point measure, hence the need to define measure-valued PDMPs. The measure-valued process representation in population dynamics is used e.g. in [START_REF] Bertoin | Random fragmentation and coagulation processes[END_REF][START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF], with fragmentation-type processes. It is a particular case of measure-valued PDMPs, with no deterministic dynamics between jumps and exponential distributions for the jump times.

After constructing measure-valued PDMPs, we define a sequence of random horizon optimal stopping problems for measure-valued PDMPs and prove that the value functions can be obtained by iterating some dynamic programming operator. We also exhibit a sequence of ε-optimal stopping times. Our approach is based on [START_REF] Gugerli | Optimal stopping of a piecewise-deterministic Markov process[END_REF] that solved the optimal stopping problem for finite-dimensional PDMPs.

When dealing with a branching population, some important characteristics of the global population, e.g. laws of large numbers for functionals of the individuals, can be obtained by simply studying a suitably weighted random tagged lineage, by means of many-to-one formulas, see e.g. [START_REF] Bertoin | Random fragmentation and coagulation processes[END_REF][START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF][START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF]. Here, we prove that this property does not hold true for the optimal stopping problem. We provide a simple counter-example of cell division where stopping a suitably chosen tagged cell and the whole population yield different value functions.

The paper is organized as follows. In Section 2, we construct measure-valued PDMPs.

In Section 3, we state and solve the optimal stopping problem for such processes. Finally in Section 4 we compare the value functions of the optimal stopping problems for the whole population and a tagged lineage.

Construction of measure-valued PDMPs

This section is dedicated to the construction of piecewise deterministic Markov processes taking values in some measure space. Our construction of measure-valued PDMPs follows the same lines as in [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF]: we first define the hybrid state space in which the process evolves, then we define the local characteristics giving the dynamics of trajectories between jumps, the jump times and the post-jump locations, and prove that a strong Markov process with such characteristics can be constructed. Finally we provide some toy example of such processes. We start by setting some notation that will be used throughout the paper.

Notation

Let d be a positive integer. We denote by B the σ-field of Borel sets on R d and B b its subset of bounded Borel sets. More generally, For any topological space E, we denote by B(E) its Borel σ-field, B(E) its set of measurable bounded real-valued functions, E its closure, and ∂E its boundary. Let M be the set of locally finite measures on (R d , B) and N = {µ ∈ M; ∀B ∈ B b , µ(B) ∈ N} be the set of locally finite point measures. Note that any µ in N can be expressed as a (possibly infinite) sum of Dirac distributions; the Dirac distribution with point mass at x, for x in R d , will be denoted by δ x .

Let C c (R d ) be the set of continuous real-valued functions with compact support on R d . For any measure µ ∈ M and function

f ∈ C c (R d ) set µf := R d f (x)µ(dx).
We endow M with the vague topology. Recall that the vague convergence for a sequence of measures (µ n ) n ⊂ M to a measure µ ∈ M is defined by

µ n f ---→ n→∞ µf, ∀f ∈ C c (R d ). We denote it µ n v ---→ n→∞ µ.
It is easily seen that that N is a closed subset of M for the vague topology.

For any real numbers a and b, a ∨ b and a ∧ bdenote the maximum and minimum respectively between a and b.

State space of measure-valued PDMPs

Let K be a finite set called modes or regimes space. For any mode v ∈ K, let E v be an open subset of M, representing the state space in mode v. The global state space is then

E = {(v, ζ) ∈ K × M | v ∈ K, ζ ∈ E v }.
We endow this set E with the σ-field E generated by sets of the form

A = {v} × A v for all Borel sets A v ∈ B(E v ).
δ x i ∈ N.
It gives a complete view of the population: all the information is contained in ζ. As time goes by, the number of individuals may increase or decrease leading to more or less terms in this Dirac sum representation but it will remain a measure in N. This representation is then easier to manipulate than a changing-dimension vector.

On the measure space M, we introduce a particular metric ρ for the vague topology in order to have a Polish space (i.e. a separable completely metrizable topological space). This property will be used in Section 2.3, for the explicit construction of the stochastic process. As shown in [START_REF] Kallenberg | Random measures, theory and applications[END_REF]Appendix], a suitable choice for ρ is constructed as follows. Let C be a countable basis of open bounded subsets in R d closed under finite unions. For all C in C, it exists a sequence (C n ) n in B b and an increasing sequence (f

C,n ) n in C c (R d ) such that f C,n ---→ n→∞ 1 C and 1 Cn f C,n 1 C . Since C is countable, the set {f C,n | C ∈ C
, n ∈ N} is also countable, and we then number f 1 , f 2 , . . . those functions. Any measure µ is completely determined by the set {µf k , k ∈ N}. Now, for all µ and µ ′ in M, we define the distance ρ by

ρ(µ, µ ′ ) := k 1 1 2 k 1 -exp(-|µf k -µ ′ f k |) .
From this metric on M, we define a metric on E related to its hybrid structure: any two points in E with different modes must be arbitrarily far away from each other. For all x = (v, ζ) and

x ′ = (v ′ , ζ ′ ) in E set ρ 0 (x, x ′ ) :=    2 π arctan(ρ(ζ, ζ ′ )) if v = v ′ ,
   v n = v for n m, ζ m+k ---→ k→∞ ζ in (E m , ρ). (2.1)
We thus denote

E = {(v, ζ) ∈ K × M | v ∈ K, ζ ∈ E v },
the closure of E for the distance ρ 0 . The following statement is then straightforward.

Lemma 2.1. The metric space (E, ρ 0 ) is a Polish space.

Construction of measure-valued PDMPs

We now introduce the three local characteristics of the PDMP specifying the dynamics of trajectories between jumps, the jump times and the post-jump locations.

• The flow Φ is defined by (

x, t) → Φ(x, t) = (v, Φ v (ζ, t)) for all x = (v, ζ
) in E and non-negative t, where the functions Φ v : M × R → M are continuous and have a semi-group property: for all s, t 0, we have

Φ v (•, t + s) = Φ v (Φ v (•, s), t).
The flow describes the deterministic trajectory of the process between jumps. Let

t * (x) = inf{t > 0, Φ v (t, ζ) ∈ ∂E v }
be the deterministic time the flow takes to reach the boundary of the domain starting from x = (v, ζ) ∈ E, with the usual convention inf ∅ = +∞. An infinite exit time t * means that the process cannot reach the boundary in finite time.

• The jump intensity λ : E → R + is a mesurable function, with a local integrability property: for all x = (v, ζ) ∈ E there exists some ε > 0 such that

ε 0 λ(v, Φ v (ζ, s))ds < ∞.
It determines the frequency of the jumps.

• The Markov kernel Q : E × B(E) → [0, 1] selects the post-jump locations. It has the following property: ∀x ∈ E, Q(x, {x}) = 0,
meaning that he process cannot have a no-move jump.

From these local characteristics, one can construct a stochastic process similarly to [10, Section 24] as follows. Let (Ω, F, P) be the canonical space for a sequence (U n ) n 1 of independent random variables with uniform distribution on [0, 1]. The sample path of an E-valued PDMP (X t (ω)) t 0 starting from a fixed initial point x = (v, ζ) ∈ E and for some ω ∈ Ω is defined iteratively. Let

F (x, t) = 1 {t<t * (x)} exp - t 0 λ(v, Φ v (ζ, s))ds (2.2) for (x, t) = (v, ζ, t) ∈ E × R + and Ψ 1 be the function from E × [0, 1] onto R + defined by Ψ 1 (x, u) = inf{t 0; F (x, t) u},
and define S 1 (ω) = T 1 (ω) = Ψ 1 (x, U 1 (ω)) the first jump time of the process. Thus F (x, •) is the survivor function of T 1 .
As the spaces M and E are Polish, one can use [18, p. 6] to obtain that {Q(x, •)} x∈E is a collection of probability measures on E, with a measurable dependence on the parameter x in E. Then there is a measurable function

Ψ 2 : E × [0, 1] → E such that the distribution of Ψ 2 (y, U ) is Q(y, •) for any random variable U with uniform distribution on [0, 1]. Hence, one can then define X t (ω) = Φ v (ζ, t), for 0 t < T 1 , X T 1 (ω) = Ψ 2 (v, Φ v (ζ, T 1 (ω))), U 2 (ω) .
Hence the trajectory X t follows the deterministic flow starting from

X 0 = x until the first jump time T 1 . At T 1 a new location Z 1 = X T 1 is drawn according to the Markov kernel Q. Namely, the law L(Z 1 | T 1 ) of Z 1 conditionally on T 1 is L(Z 1 | T 1 ) = Q(Φ(x, S 1 ), •), (2.3)
The process now restarts from X T 1 following the same steps. Define

S 2 (ω) = Ψ 1 (X T 1 (ω), U 3 (ω)), T 2 (ω) = S 1 (ω) + S 2 (ω),
and set

X t (ω) = Φ(X T 1 (ω), t -T 1 (ω)), for T 1 (ω) t < T 2 (ω), X T 2 (ω) = Ψ 2 (v, Φ(X T 1 (ω), S 2 (ω))), U 4 (ω) ,
and so on. In order to avoid explosion issues, the following assumption will hold throughout the paper.

Assumption 2.1. For all (x, t) in E × R + , E x [ ∞ n=1 1 (Tn<t )] < ∞.
Hence the trajectories of (X t ) are well defined for all t 0. The positive random variables S 1 , S 2 , . . . are the times between two consecutive jumps or inter-jump times. For notational convenience, we set T 0 = S 0 = 0 and T n := n i=1 S n the n th jump time of the process. Note that we have T 1 = S 1 and S n+1 = T n+1 -T n . The sequence (Z n ) n with Z n = X Tn describes the post jump locations of the process. By construction, all the randomness of the continuous-time process (X t ) t 0 is contained in the discrete-time process (Z n , S n ) n 0 .

Structure of stopping times and Markov property

The aim of this section is to prove that the special structure of stopping times for finite dimensional PDMPs given in [START_REF] Davis | Markov models and optimization[END_REF]Theorem A2.3] still holds in the measure space context. This yields the Markov property and will be important for the study of the optimal stopping problem in Section 3.

Let D E [0, +∞) be the set of unctions on R + with values in E that are rightcontinuous with left limits. Denote by Xt the coordinate function Xt (f ) = f (t) for f ∈ D E [0, +∞). Let (F 0 t ) t 0 denote the natural filtration of (x t ) and F 0 = ∨ t 0 F 0 t . Under Assumption 2.1, for each starting point x ∈ E, the construction in the previous section defines a measurable mapping Ψ x from Ω onto D E [0, +∞) such that Xt (Ψ x (ω)) = X t (ω). Let P x denote the image measure of P by Ψ x . This defines a family of measure (P x ) x∈E on D E [0, +∞). In the sequel, we identify Xt and X t . For any probability measure ν on E define the measure P ν on (D E [0, +∞), F 0 ) by

P ν (• • • ) = E P x (•)ν(dx). Now let F ν
t be the completion of F 0 t with all P ν -null sets of F 0 and define F t = ∩ ν∈P(E) F ν t , where P(E) is the set of probability measures on E.

We can now state a crucial result on the structure of stopping-times for our process.

τ = ∞ n=1 R n-1 ∧ S n .
Proof. The proof follows the same lines as in [START_REF] De Saporta | Numerical methods for simulation and optimization of piecewise deterministic Markov processes: application to reliability[END_REF]Section 1.7]. It is based on Theorem A2.3 in [START_REF] Davis | Markov models and optimization[END_REF] that is valid for any right-continuous piecewise constant process taking values in a Borel set and the one-to-one correspondence between our process (X t ) and the right-continuous piecewise constant E × N-valued process (η t ) defined by

η t = (X 0 , 0), t < T 1 , η t = (X Tn , n), T n t < T n+1 ,
see [10, (25.1)].

Theorem 2.2. The process (X t ) t 0 on (D E [0, +∞), F, (F t ), (P x ) x∈E ) is a strong Markov process.

Proof. The strong Markov property is proved in the same way that in [10, Section 25], using Theorem 2.1 and the one-to-one correspondence with the piecewise constant process (η t ) defined above.

Toy example

We now develop Example 2.1 into a more generic model for some cells population. Typically, cells grow and divide into two daughter-cells that start growing and then divide and so on.

Single cell model

We first define a single-cell model that follows standard final dimensional PDMP dynamics, where one randomly selects a single daughter cell at each division. This model was further studied in [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF]. We consider a model with a single mode, hence the state space is simply R + . For ξ in R + and t > 0, the size at time t for a cell with initial size ξ is given by φ(ξ, t) = ξ exp(rt), where r > 0 is the common growth rate for all the cells. The jump intensity giving the division dates is given by l(ξ) = ξ α , for some positive number α. This is simplified model consistent with the statistical evidence that division is triggered by the cell size rather than its age, see [START_REF] Robert | Division in escherichia coliis triggered by a size-sensing rather than a timing mechanism[END_REF]. The jump kernel is simply a division by two Q(x, A) = δ x/2 (A) for any Borel subset A of R + .

Population model

We now consider the previous cell growth-division dynamics but for the whole population instead of a single randomly selected cell. Again we consider a single mode so that he space E is simply N.

Let ζ = n i=1 δ x i be in N be an initial state of n cells with respective sizes (x i ) 1 i n . Each cell grows following the previous dynamics, so that globally the flow is

Φ : N × R + → N (ζ, t) → n i=1 δ φ(x i ,t) .
The first jump time corresponds to the first split ie the minimum between n exponentially distributed random variables. Thus, the jump intensity is:

λ : N → R ζ → n i=1 l(x i ).
Given that the j th cell is the one which split, the post-jump location is

Z 1 = n i=1 i =j δ x i + 2δ x j /2 ,
where celle number j was removed and two new cells are added with half the size of cell j. For notational convenience, we set

ζ (j) := ζ -δ x j + 2δ x j /2 .
Thus, for A in E, the Markov jump kernel of the PDMP is given by

Q(ζ, A) = n j=1 l(x j ) n i=1 l(x i ) 1 A (ζ (j) ).

Time-augmented population model

Sometimes, it is convenient to add time in the state variable of a PDMP, for instance when one wants to use time-dependent jump intensity, or trigger some jump when a certain lapse of time has passed, or to study control problems with time-dependent reward functions. In [10, Section 31] Davis proves that for finite dimensional PDMPs, the time-augmented process is still a PDMP. The same property holds in the framework of our toy example. We consider the time-augmented process X t : defined as the PDMP on N starting from the initial state ζ = n i=1 δ (x i ,0) . The flow is now given by

Φ : N × R + → N ( ζ, t) → n i=1 δ (φ(x i ,t),u+t) ,
for any ζ in N of the form ζ = n i=1 δ (x i ,u) . The jump intensity is

λ : N → R ζ = n i=1 δ (x i ,u) → n i=1 l(x i ),
and the jump kernel is simply

Q( n i=1 δ (x i ,u) , A) = n j=1 l(x j ) n i=1 l(x i ) 1 A ( n i=1 δ (x i ,u) -δ (x j ,u) + 2δ (x j /2,u) ).
We will further study this toy example in Section 4 to prove that controlling the whole population is not equivalent to controlling a suitably chosen random lineage.

Optimal stopping problem

We now turn to the main aim of this paper: defining and solving the optimal stopping problem for measure-valued PDMPs. Roughly speaking, one wants to stop the process at the best time in order to maximize some reward depending on the state of the process when stopped. More precisely, let (X t ) t 0 be an E-valued PDMP on (D E [0, +∞), F, (F t ) t 0 , (P x ) x∈E ) and g ∈ B(E) be some non-negative reward function. Denote by M the set of stopping-times with respect to the filtration (F t ) t 0 and for all positive integer N , let M N be the set of stopping-times bounded by the N th time jump

T N of the PDMP M N = {τ ∈ M; τ T N }. For all x ∈ E, set V(x) := sup τ ∈M N E[g(X τ ) | X 0 = x]. (3.1)
Thus V(x) is the best possible (average) performance when stopping a PDMP starting from X 0 = x before its N -th jump. Function V is called the value function of the optimal stopping problem. Solving an optimal stopping problem consists in characterizing the value function as the solution of some recursive equations called dynamic programming equations and exhibiting a family of ε-optimal stopping times τ ε ∈ M N such that

V(x) -ε E[g(X τε ) | X 0 = x] V(x).
In this section, we first define some suitable dynamic programming operators and a family of stopping times. Then we prove that the value functions can be constructed by iteration of the dynamic programming operators and that the stopping times are ε-optimal. This section is inspired from the study of the optimal stopping problem for finite (fixed) dimension PDMPs derived in [START_REF] Gugerli | Optimal stopping of a piecewise-deterministic Markov process[END_REF].

Dynamic programming operators

We start with some additional notation and assumption. For w in B(E), x in E and l a measurable real-valued function on E, we denote in short

lQw(x) := l(x) × Qw(x) = l(x) E w(y)Q(x, dy).
The following assumption is made for simplicity reasons. It is satisfied in most real-life examples for instance when monitoring a population until some finite horizon time.

Assumption 3.1. The exit time t * is in B(E).

We now define some operators on B(E). Let H and I be the operators from B(E) onto B(E × R + ) defined for all w in B(E), x in E and t ∈ R + , by

Hw(x, t) = w (Φ(x, t ∧ t * (x))) e -Λ(x,t∧t * (x)) , Iw(x, t) = t∧t * (x) 0 λQw(Φ(x, s)) e -Λ(x,s) ds, where Λ(x, t) = t 0 λ(Φ(x, s))ds.
We also introduce operator K from B(E) onto B(E), defined for all w in B(E) and x in E as

Kw(x) = t * (x) 0 λQw(Φ(x, s)) e -Λ(x,s) ds + Qw(Φ(x, t * (x)))e -Λ(x,t * (x)) .
It is straightforward to see that these operators can be expressed as expectations involving the embedded Markov chain (Z n , S n ) defined in Section 2.3. Proposition 3.1. For all w in B(E), x in E and t 0 one has:

Hw(x, t) = E x w(X t∧t * (x) ) 1 S 1 >t∧t * (x) = w (Φ(x, t ∧ t * (x))) P x (S 1 > t ∧ t * (x)) , Iw(x, t) = E x w(Z 1 ) 1 S 1 t∧t * (x) , Kw(x) = E x [w(Z 1 )].
Finally we denote by J and L the operators from B(E) onto B(E × R + ) and B(E) respectively defined for all w in B(E), x in E and t ∈ R + by J(w, g)(x, t) = Hg(x, t) + Iw(x, t),

L(w, g)(x) = sup t 0 {J(w, g)(x, t)} ∨ Kw(x),
where g is the reward function of the optimal stopping problem. Roughly speaking, operator L represents the best compromise between stopping at the best location along the deterministic trajectory following the flow (sup t 0 J part) or waiting until the next jump (K part).

Family of stopping-times

Now we introduce a family of random variables and prove they are stopping times. They will be candidates ε-optimal stopping times for our optimal stopping problem. For all x ∈ E, n ∈ N and ε > 0 set

r n,ε (x) = t * (x) if KV n (x) > sup t>0 J(V n , g)(x, t), inf{s 0 ; J(V n , g)(x, s) sup t>0 J(V n , g)(x, t) -ε} otherwise.
For all ε > 0 and n 2 also set

R 1,ε = r 0,ε (Z 0 ), R ε n,0 = r n-1,ε/2 (Z 0 ), R ε n,k = r n-k-1,ε/2 k (Z k )1 (R ε n,k-1 S k ) , 1 k n -1. Then R ε n,k is clearly F T k -measurable for 0 k n -1. Finaly, define S 1,ε := r 0,ε (Z 0 ) ∧ T 1 = R ε 1,0 ∧ S 1 ,
and by iteration,

S n,ε := R ε n,0 T 1 > R ε n,0 , T 1 + θ(T 1 )S n-1,ε/2 T 1 R ε n,0 ,
where θ(t) is the shift operator with lag t on D E [0, +∞), namely for

f ∈ D E [0, +∞), θ(t)f (•) = f (t + •).
In order to prove that the S n,ε are stopping-times in M n , we first study the effect of the shift operator on R ε n,k . Lemma 3.1. For all ε > 0 and n 2 and 1 k n -1, on the set (T 1 S n,2ε ), one has

R 2ε n,k = θ(T 1 )(R ε n-1,k-1 ). Proof. For n = 2, by definition, one has R ε 1,0 = r 0,ε (Z 0 ) hence θ(T 1 )(R ε 1,0 ) = r 0,ε (Z 1 ) and R 2ε 2,1 = r 0,ε (Z 1 ) on (T 1 S 2,2ε ) = (T 1 R ε 2,0
). Hence the result holds. For n 3 we prove the result by induction on k using similar arguments and the fact that (T 1 S n,ε ) = (R 2ε n,0 T 1 ).

We now prove that S n,ε is a stopping-time using the characterization of Theorem 2.1. Lemma 3.2. For all ε > 0 and n ∈ N * one has

S n,ε = n k=1 R ε n,k-1 ∧ S k .
In particular, S n,ε is a stopping time and S n,ε T n .

Proof. We proceed by induction on n. For n = 1, by definition

S 1,ε = r 0,ε (Z 0 ) ∧ T 1 = R ε 1,k-1 ∧ S 1
and the result is true. Suppose the result holds for n -1. From the definition, on the event (r

n-1,ε/2 (Z 0 ) = R ε n,0 < S 1 ), one has S n,ε = R ε n,0 and R ε n,1 = 0 < S 2 , R ε n,2 = 0 < S 3
, and so on, so that

n k=1 R ε n,k-1 ∧ S k = R ε n,0 and the result is valid. On the event (R ε n,0 S 1 ), by definition S n,ε = T 1 + θ(T 1 )S n-1,ε/2
. Now, the induction hypothesis and the previous lemma yield

T 1 + θ(T 1 )S n-1,ε/2 = S 1 + n-1 k=1 R ε n,k ∧ S k+1 = n k=1 R ε n,k-1 ∧ S k , on (R ε n,0 S 1 )
. Hence the result.

Characterization of the value function

We can now propose an iterative construction of value functions by iterating operator L. For all x ∈ E, set

V 0 (x) = g(x), V n (x) = L(V n-1 , g)(x) for all n 1. (3.2)
Clearly the functions V n are in B(E). We can now state and prove our main result, namely that the value function V of the optimal stopping problem (4.2) equals the N -th iterate V N and that S N,ε is an ε-optimal stopping time.

Theorem 3.1. Let x be in E and n ∈ N. Then one has, for all ε > 0, S n,ε is a stopping-time in M n and

V n (x) = sup S∈Mn E x [g(X S )], (3.3)
E x [g(X Sn,ε )] V n (x) -ε. (3.4)
Proof. By an induction argument we will prove Equation (3.4) and the following inequality

∀S ∈ M, ∀n ∈ N, E x [g(X S∧Tn )] V n (x). (3.5)
These two equations imply (3.3). Indeed, for all S in M ∞ , we have S ∧ T n ∈ M n then for all S in M ∞ , n in N and x in E, (3.4) and (3.5) yield

V n (x) -ε E x [g(X Sn,ε )] sup S∈Mn E x [g(X S )] = sup S∈Mn E x [g(X S∧Tn )] V n (x),
which is valid for all positive ε, hence the result.

It remains now to prove (3.5) and (3.4) by induction.

The case n = 1 is based on Theorem 2.1. Indeed, we obtain E x [g(X S∧T 1 )] = E x [g(X R 0 ∧T 1 )], for some F 0 -measurable random variable R 0 . From this, we deduce

E x [g(X S∧T 1 )] = E x [g(X R 0 )1 T 1 >R 0 ] + E x [g(Z 1 )1 T 1 R 0 ] = Hg(x, R 0 ) + Ig(x, R 0 ) = J(g, g)(x, R 0 ) sup t 0 {J(g, g)(x, t)} V 1 (x);
which proves (3.5) for n = 1. To prove (3.4), we distinguish two cases:

• if Kg(x) > sup t 0 J(g, g)(x, t), then V 1 (x) = Kg(x) = E x [g(Z 1
)] and also

S 1,ε = t * (x) ∧ T 1 = T 1 . Thus E x [g(X S 1,ε )] = E x [g(Z 1 )] = V 1 (x) V 1 (x) -ε. • otherwise, V 1 (x) = sup t 0 J(g, g)(x, t) and, by definition of J, E x [g(X S 1,ε )] = E x [g(X r 0,ε (x)∧T 1 )] = E x [g(X r 0,ε (x) )1 T 1 r 0,ε (x) ] + E x [g(Z 1 )1 T 1 r 0,ε (x) ]
= J(g, g)(x, r 0,ε (x)).

We then deduce, by definition of r 0,ε :

E x [g(X S 1,ε )] sup t 0 {J(g, g)(x, t)} -ε = V 1 (µ) -ε.
This completes the proof of (3.4) for n = 1. Now, let N 1, suppose that (3.5) and (3.4) hold true for all n N and prove that they hold for n = N + 1.

Begin by proving (3.5). Again, this is based on Theorem 2.1. As S is a stopping-time, it can be decomposed as

S = ∞ n=1 R n-1 ∧ S n .
where (R n ) n∈N are non-negative (F Tn ) n∈N -adapted random variables. In particular, on (S T 1 ) one has

S = S 1 + ∞ n=2 R n-1 ∧ S n = T 1 + θ(T 1 )(S ′ ),
for some stopping-time S ′ . Thus one has

E x [g(X S∧T n+1 )] = E x [g(X S∧T 1 )1 S<T 1 ] + E x [g(X S∧T n+1 )1 S T 1 ] = E x [g(X R 0 ∧T 1 )1 R 0 <T 1 ] + E x [g(X S∧T n+1 )1 R 0 T 1 ] = E x [g(X R 0 )1 R 0 <T 1 ] + E x [E x [g(X (T 1 +θ(T 1 )(S ′ ))∧T n+1 )|F 1 ]1 R 0 T 1 ] = E x [g(X R 0 )1 R 0 <T 1 ] + E x [E Z 1 [g(X S ′ ∧Tn )]1 R 0 T 1 ] (3.6) E x [g(X R 0 )1 R 0 <T 1 ] + E x [V n (Z 1 )1 R 0 T 1 ] (3.7) sup r 0 E x [g(X r )1 r<T 1 ] + E x [V n (Z 1 )1 r T 1 ] sup r 0 J(V n , g)(x, r) ∨ KV n (x) = V n+1 (x).
Line (3.6) is obtained thanks to the strong Markov property. The induction hypothesis is applied in line (3.7). This achieves the induction for (3.5).

Let us now prove (3.4). One has

E x [g(X S n+1,2ε )] = E x [g(X S n+1,2ε )1 R ε n,0 <T 1 ] + E x [g(X S n+1,2ε )1 R ε n,0 T 1 ] = E x [g(X R ε n,0 )1 R ε n,0 <T 1 ] + E x [g(X T 1 +θ(T 1 )Sn,ε )1 R ε n,0 T 1 ] = E x [g(X R ε n,0 )1 R ε n,0 <T 1 ] + E x [E Z 1 [g(X Sn,ε )]1 R ε n,0 T 1 ] E x [g(X R ε n,0 )1 R ε n,0 <T 1 ] + E x [V n (Z 1 )1 R ε n,0 T 1 ] -ε × P x (R ε n,0 T 1 ) E x [g(X R ε n,0 )1 R ε n,0 <T 1 ] + E x [V n (Z 1 )1 R ε n,0 T 1 ] -ε. From the definition of R n,ε one readily obtains • if KV n (x) > sup t 0 J(V n , g)(x, t), then V n+1 (x) = KV n (x) = E x [V n (Z 1 )] and R ε n,0 = t * (x). So the previous inequality becomes E x [g(X S n+1,2ε )] E x [g(X t * (x) )1 t * (x)<T 1 ] + E x [V n (Z 1 )1 T 1 t * (x) ] -ε = KV n (x) -ε = V n+1 (x) -ε.
• Otherwise, V n+1 (x) = sup t 0 J(V n , g)(x, t), and the previous inequality becomes

E x [g(X S n+1,2ε )] J(V n , g)(x, r 0,ε (x)) -ε sup t 0 J(V n , g)(x, t) -ε -ε = V n+1 (x) -2ε.
Thanks to these two cases, we prove (3.4) and we end the proof.

Comparaison between the tagged cell and measure-valued process

In this section, we investigate wether stopping a single well chosen individual is equivalent to stopping the whole population. When dealing with a branching population, some important characteristics of the global population, e.g. laws of large numbers for functionals of the individuals, can be obtained by simply studying a suitably weighted random tagged lineage, by means of many-to-one formulas, see e.g. [START_REF] Bertoin | Random fragmentation and coagulation processes[END_REF][START_REF] Cloez | Limit theorems for some branching measure-valued processes[END_REF][START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF][START_REF] Guyon | Limit theorems for bifurcating Markov chains. Application to the detection of cellular aging[END_REF]. We prove that this property does not hold true for the optimal stopping problem. We use the simple toy example of cell division from Section 2.5 to show that stopping a suitably chosen tagged cell and the whole population yield different value functions.

Tagged cell and many-to-one formula

Let us rapidly describe the many-to-one formula and the definition of the tagged as presented in [START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF]. Heuristically, picking a branch uniformly at random along the genealogical tree describe the path of a tagged cell whose the behaviour is similar to the one of an individual picked at time t. More precisely, for any measurable positive function h and any t 0 we have the following formula, commonly called many-to-one formula

E x u∈U X t (h)e -rt = n i=1 E x i h(χ t ) χ t × x i , ( 4.1)
where x = µ = n i=1 δ (x i ,0) , (X t ) t 0 is the measure-valued PDMP of Section 2.5.2, and (χ t ) t 0 is the tagged process, representing the evolution of the size of the tagged cell. Its dynamics is as a real-valued PDMP described in Section 2.5.1 but whose parameters are different from those of (X t ) t 0 ; see [START_REF] Bertoin | Random fragmentation and coagulation processes[END_REF][START_REF] Cloez | Limit theorems for some branching measure-valued processes[END_REF][START_REF] Doumic | Statistical estimation of a growth-fragmentation model observed on a genealogical tree[END_REF] for details.

To compare the value functions of the tagged cell and the measure-valued process, we will impose the form of our reward functions based on Equation 4.1. For the tagged cell, we choose a bounded nonnegative function f continuous along the flow. For the process (X t ) t , the reward function g has the form g :

n i=1 δ (x i ,t) ∈ E → n i=1 f (x i )e -rt , ( 4.2)
where we used the time-augmented process defined in Section 2.5.3 to take the time dependence into account.

Comparison of the value functions

Numerically computing the value function V N is very demanding, as for each iteration of the dynamic programming operators, one needs to compute the functions on the whole state space E, at least at first sight. Actually, when dealing with the optimal stopping problem with horizon N , the dynamic programming recursion on functions V n can be rewritten as a recursion on the random variables V k (Z N -k ). Still, the recursion is numerically intractable as one needs to compute conditional expectations at each step.

In order to avoid such intricacies, we simply consider the optimal stopping problem with horizon N = 1 jump.

We thus have to compute

V 0 (Z 1 ) = g(Z 1 ), V 1 (Z 0 ) = sup t 0 {E Z 0 [V 0 (Z 1 )1 T 1 t ] + g • Φ(Z 0 , t) × P Z 0 (T 1 > t)} ∨ E Z 0 [V 0 (Z 1 )].
The only quantities to look at are V 0 (Z 1 ) = g(Z 1 ) and V 1 (Z 0 ) = V 1 (µ). Moreover, this last expression is deterministic if we choose a deterministic Z 0 . More specifically, one has

V 0 (Z 1 ) = e -rT 1   i =I 1
(x i e rT 1γ)1 x i <γe -rT 1 + (x I 1 e rT 1 -2γ)1 x I 1 <2γe -rT 1 + (n + 1)γ   , (4.3) and

V 1 (Z 0 ) = sup t>0    E x   e -rt   i =I 1
(x i e rT 1γ)1 x i <γe -rT 1 + (x I 1 e rT 1 -2γ)1 x I 1 <2γe -rT 1 + (n + 1)γ

  1 T 1 t   + e -rt i (x i e rt -γ)1 x i <γe -rt + nγ exp -(x α i ) e rtα -1 αr ∨ E x [V 0 (Z 1 )], (4.4
) where the first jump time T 1 is distributed as (2.2) and the random variable I 1 is the rank of the split cell. Its distribution is given by P µ (I 1 = j) = . We can numerically simulate V 0 (Z 1 ). For V 1 (µ), we fix a (large enough) maximum time and discretize the time interval in order to compute (an approximation of) the supremum. In the same way, the value function for the tagged cell is given by: V 1 (x) = sup t 0 (xe rtγ)1 xe rt <γ + γ exp -

x α αr (e rαt -1)

+ E x (xe rT 1 /2 -γ)1 xe rT 1 <2γ + γ 1 T 1 t ∨ E x [V 0 (Z 1 )] .
For the numerical results, we use the following parameters: α = 1, r = 2, γ = 1. For the discretization step, we evaluate the maximum value of T 1 by the Monte-Carlo simulations and we divide it by nbpt which corresponds to the number of discretization points for the evaluation of the supremum. Here, nbpt = 10000 and the Monte Carlo sample size is N = 100000. We obtain the following approximation V 1 (3) = 1.0018 ± 2.54 × 10 -4 V 1 (δ 3 ) = 1.3447 ± 3.6034 × 10 -4 .

(4.5)

Hence, even on this simple toy example where the cost functions g and h match and with a short one-step horizon, the optimal performance for the global population and the tagged cell differ. It will be the object of future work to design specific numerical approximations of the value function of the global measure-valued population.

  For instance, if one can consider the temporal evolution of a cell population characterized by state variables such as age, size, growth, maturity, protein content... Then the quantity ζ ∈ N represents such state variables and the mode can be the experiment conditions. In the following example, we consider a simple case with only one mode and the variable state corresponds to the cell size.

	Example 2.1. Consider a population of 3 individuals at a given time. Their sizes are
	denoted by (x i ) 1 i 3 . This population is identified with the measure
	3
	ζ =
	i=1

  (x, x ′ ) equals 1 if and only if x and x ′ have different modes. With this metric, a sequence (x n ) n = (v n , ζ n ) n converges to x = (v, ζ) in (E, ρ 0 ) if and only if it exists some m in N such that:

	Thus ρ 0
	1	otherwise.

* This work was partially supported by Région Languedoc-Roussillon and FEDER under grant Chercheur(se)s d'Avenir, project PROMMECE


Jalal Possik
email: jpossik@yorku.cacharlesyaacoub@usek.edu.lb

Charles Yaacoub

Simon Gorecki
email: simon.gorecki@u-bordeaux.fr

Gregory Zacharewicz Gretha
email: gregory.zacharewicz@mines-ales.fr

Andrea D ' Ambrogio

A MODEL BASED SYSTEMS ENGINEERING APPROACH TO AUTOMATED MUSIC ARRANGEMENT

Keywords: model-based systems engineering (MBSE), automated music arranger, music generation, BPMN, MIDI parsing

MBSE (Model Based Systems Engineering

) is considered a valuable and effective approach not only in engineering domains, although it appears to bring a potential assistance to artistic visions and initiatives. The cross-domain capitalization is the key to further innovation and advances in technology. This paper describes an MBSE approach, applied to the development of a complex multi-component arrangement system, to generate an automated music arrangement for any melody, with a harmonically correct result. The user introduces the melody in the form of a MIDI file. The system analyzes the melody, detects the scale, then follows algorithms based on thorough theoretical and harmonic studies, to finally generate the arrangement based on the selected music genre (classical, jazz, pop, etc.). Compared to other existing systems in the market, the developed system largely reduces the dissonance between the melody and the arrangement, which makes the final arrangement a good accompaniment to the melody.

INTRODUCTION

The digitalization of concepts, ideas, and propositions introduced to satisfy needs that may arise in any social or business domains, contributes to enhance creativity and provide added value. When moving to a digital world, a major problem frequently identified is the transformation of human-related vision and senses into digitalized systems. The contribution in [START_REF] Zacharewicz | Model-Based Approaches for Interoperability of next Generation Enterprise[END_REF] clearly identifies the opportunities provided by modeling as an effective approach to achieve such an objective. A model may represent a system, a structure, a phenomenon, or a concept. Researchers and engineers are increasingly using the model-based systems engineering (MBSE) approach to describe and control complex systems. Moreover, MBSE and models simulation techniques are widely applied to examine "what if" scenarios that help managers in decision making. In this paper, the MBSE approach is employed as a guide to the development and implementation of a complex automated music arranger system that organizes the music according to the needs or requirements of the composer, the group of performers, or the singer.

To compose a piece of music, the melody is created first, then the musical arrangement is produced. Depending on the style, the music arranger makes the arrangement, i.e., the arrangement of classical music is different from the arrangement of jazz music. The instruments used, the music rules, and the rhythm differ. A musical arrangement can be played by a band or an orchestra, although, it can also be obtained by using synthesizers or Virtual Studio Technology (VSTi) instruments. Musical arrangements always requires a human intervention; no machine, to our knowledge, has been able to generate a correct chord progression to accompany a certain melody. The arrangement process of a three-minute melody takes hours, and based on the chosen style, it might take days sometimes. The goal of the project described in this paper is to simplify the work of music arrangers by developing a multicomponent system that parses the MIDI melody file and converts it to a fully arranged music within seconds.

In summary, this project aims at fully automating the music arrangement process using a model-based approach, while respecting the rules of music harmony. After providing the melody to the developed system, the user chooses a set of parameters (musical genre, tempo, etc.). Then, the system analyzes the melody, associates different chord progressions with it, and finally performs the simulation process to select the appropriate chord progression that accompanies the melody.

STATE OF THE ART

Research efforts have been made to combine artistic activities with the current digital world revolution. [START_REF] Bonnin | Automated Generation of Music Playlists: Survey and Experiments[END_REF] presented a comparative evaluation of playlist generation methods based on historical data. [START_REF] Morris | Music Generation Using Cellular Models[END_REF] introduced a mapping tool for the automated creation of music using cellular models. [START_REF] Whorley | Music Generation from Statistical Models of Harmony[END_REF] investigates an improved sampling technique for music generation from statistical models of harmony. Nevertheless, most of the music generation research involves both the melody and the arrangement generation. When the melody is imposed by the user, the generation of arrangement becomes more difficult and complicated requiring deep theoretical and harmonic studies of the melody to detect the appropriate chords (for each beat/bar) on which an arrangement is built, based on the selected music genre (classical, jazz, pop, etc.).

Most of the existing arranger tools are non-automatic arrangers. In this kind of software, there is no scale detection, and subsequently, no automatic generation of chords. The user provides the chords he sees fit, and the arrangement will be generated according to these chords. These software tools are based only on user chords and arrangement styles provided by the arranger, without any theoretical analysis performed by the machine. There are, to our knowledge, only two tools for automatic arrangement, in addition to machine learning-based tools, but not based on theoretical and harmonic studies. This leads to wrong scale detection as well as dissonant chords with the melody. This section presents some of the existing arranger software, highlighting their ineffectiveness from a professional's point of view, which motivates our idea of developing an automatic arranger. For non-automatic arrangement, several software tools exist. For instance, "Easy Music Composer" [START_REF] Biglobe | Easy Music Composer (version 9.97[END_REF] and "Chord Composer" (kvraudio 2021). In such tools, arrangement is generally obtained after three steps: loading the melody in MIDI file format, choosing the arrangement style, and manually providing the suitable chords for the loaded melody. Such tools might be convenient whenever the user's own chords are desired. On the other hand, "Onyx Arranger" (Jasmine Music Technology 2004) and "e-Xpressor" [START_REF] Maruta | E-Xpressor (version 1.31[END_REF] allow the automated music arrangement after loading the melody as a MIDI file and choosing the desired style. With both tools, in many cases the arranged music has a lot of harmonic errors, which most of the time cause dissonance. Moreover, it is frequently observed that silence is suddenly introduced, i.e. the arrangement suddenly disappears. Furthermore, frequent chord repetition is also observed. The "Onyx Arranger" also suffers from a poor scale detection, which leads to complete dissonance throughout the melody. Recently, "AWS DeepComposer" (Amazon 2020) has emerged as a new arranger tool based on machine learning approaches. The user provides the tool with a simple melody which is sent to a cloud service provided by Amazon Web Services, where generative artificial intelligence is exploited not only to arrange the melody, but also to complete the melody with machine-generated music. The synthesized music includes guitar pieces, bass, synths and drums. Generative models are provided for four different genres: rock, pop, jazz, and classical, in addition to the user's own models that can be built. The main disadvantage of DeepComposer is that it highly relies on training length and datasets; larger datasets and extensive trainings often provide better results, at the expense of increased computational resources and costs. In addition, theoretical correctness is highly dependent on training data that require hundreds of musical pieces for each genre, which might not be convenient for users.

The current project opposes the artistic view on the one hand, and the technical view on the other hand. Some research efforts have been produced to bridge the gap between both different worlds, for instance, model-based reconciliation can be considered one of the important solutions that can be used (D'Ambrogio and Bocciarelli 2007). Model driven engineering methods have emerged to facilitate the trans-domain interoperability and support productive modeling. The first and most used one is the Model Driven Architecture (MDA) proposed by the Object Management Group (OMG 2013). MDA is a system design approach that introduces a significant degree of automation in the development of software systems. It provides a set of guidelines for structuring and transforming specifications, which are expressed as models. It is considered as an important solution in diverse research domains [START_REF] Gianni | Modeling and Simulation-Based Systems Engineering Handbook[END_REF], especially for the prediction, analysis, and description of complex systems [START_REF] Possik | A BPMN/HLA-Based Methodology for Collaborative Distributed DES[END_REF]). Nevertheless, model-based approaches have not been used in the music domain yet. In this paper, we present a model-based approach to drive the development and implementation of a complex multicomponent music arrangement system. The Business Process Model and Notation (BPMN) is used to define the existing system components and their interaction. Papyrus, an integrated environment for editing UML (Unified Modeling Language) models and executing UML-based processes, is used to accomplish the orchestration between the developed components.

MATERIALS AND METHODS

MIDI (Musical Instrument Digital Interface

) was developed in 1983 as a software and hardware design that enables the communication between synthesizers and computers. The MIDI communication protocol allows the data exchange between electronic musical instruments. The MIDI digital interfaces are designed to transmit control information, called messages or events, between devices. Unlike audio recordings, which capture sound waves, MIDI recordings capture MIDI messages, which in turn broadcast the sound to an external device (sound generator), such as a synthesizer. Using MIDI format, it is possible to apply nondestructive transformations, such as transposition or tempo modifications. Such transformations only modify the MIDI messages and not the sound quality of the recording. Using MIDI files allows the users to change many factors that look hard and sometimes impossible with audio files. In our project, the melody is supplied in MIDI format, as in other existing software solutions on the market. The first component developed for the music arranger system is a MIDI parser.

𝐵𝑃𝑀

MIDI Parsing

The MIDI file is a succession of blocks called chunks. Each chunk follows the pattern of the chunk header. A MIDI file contains several chunks that can be of different sizes. Figure 1 is an example of the music notation and the hexadecimal code resulted from the recorded melody of the popular nursery rhyme "twinkle twinkle little star" in MIDI format. The hexadecimal code starts with the MIDI track header (MThd), 0x4d 0x54 0x68 0x64, that marks the beginning of a MIDI file. The length in bytes of the rest of the chunk header, the MIDI type (0,1, or 2), the number of tracks, and the number of pulses per quarter note (PPQN) are found in the rest of the hexadecimal code of the first row (see Figure 1). 0x4d 0x54 0x72 0x6b of the MIDI track (MTrk) marks the beginning of the track. The tempo hexadecimal code starts with the byte 0xff that indicates the beginning of a meta message; the second byte (0x51) represents the tempo meta message. The third byte is the size of the tempo data (0x03). Hence, the following 3 bytes (0x0b71b0) represent the number of microseconds per beat. It is a 24-bit unsigned, big-endian integer that is equal to 750,000 microseconds per beat when converted to decimal. In all applications and synthesizers, the tempo is shown in Beats Per Minute (BPM). There are 60,000,000 microseconds in a minute. The formula used to calculate the tempo in BPM is: 60,000,000 = 𝑚𝑖𝑐𝑟𝑜𝑠𝑒𝑐𝑜𝑛𝑑𝑠 𝑏𝑦 𝑏𝑒𝑎𝑡 Therefore, the tempo used in the MIDI melody of Figure 1 is equivalent to: 60,000,000 750,000 = 80 𝐵𝑃𝑀 If no tempo meta messages exist in the MIDI file, 500,000 microseconds per beat (120 BPM) are considered as the default tempo of the song. The hexadecimal code 0xff marks the beginning of a meta message. The hexadecimal 0x58 marks the beginning of the time signature related section. This byte is followed by 5 bytes. The second byte represents the numerator of the time signature. In this case, the numerator is 4. The third byte represents the denominator as a negative power of 2; i.e., 2 signifies a quarter note, 3 signifies an eighth note, etc.). In this case (0x02), it represents a quarter note. The code 0xff 0x59 marks the beginning of the key signature section. This code is followed by 3 bytes. The second byte specifies the number of flats or sharps that identify the key signature. In the case of Figure 1, it is equal to 0x00, indicating that no flats or sharps exist in the key signature. The third byte specifies a major (0x00) or minor (0x01) key.

After the hexadecimal meta messages, the synthesizer starts reading the notes. 0x90 means that a note is turned on. 0x80 means that a note is turned off. Both hexadecimal bytes (0x90 and 0x80) are followed by two bytes. The first byte specifies the note and the second one specifies the velocity of the note, i.e., 0x90 0x3c 0x4c means that the note number 60 or C4 (0x3c) is turned on with a velocity equal to 76 (0x4c). This code is followed by 0x80 0x3c 0x00, which means that the same note (0x3c) is turned off with a velocity 0. Turned on notes should be turned off. Check Figure 1 for the detailed description of the hexadecimal bytes related to the notation process. The hexadecimal sequence 0xff 0x2f 0x00 marks the end of the track.

BPMN is a business process-modeling standard that offers a graphical notation based on a flowcharting technique. BPMN represents the end-to-end flow of a process. The Business Process Management Initiative (BPMI) developed the Business Process Modeling standard. In 2005, this group merged with the OMG (Object Management Group). In 2011, OMG released BPMN version 2.0 and changed the name of the standard to Business Process Model and Notation. The modeling language became more detailed by using a richer set of symbols and notations for business process diagrams. The main goal of BPMN is to deliver a standard notation easily readable by non-expert users. In the present work, we used BPMN to specify the proposed methodology and simplify the understanding of the integration and collaboration between discrete event simulators.

Figure 1: The music notation and the MIDI hexadecimal code of the melody "twinkle twinkle little star".

The BPMN model of Figure 2 describes the MIDI parsing system developed to extract the key information required to build the automated music arranger framework. This system is composed of a MIDI parser and a header parser (each of them having a different role), in addition to Papyrus system that orchestrates the components. The parsing system loads the MIDI file from the Papyrus engine (F2. A), then reads the first four bytes to verify that the file loaded is a MIDI file. A MIDI file should start with the following sequence of hexadecimal bytes: 0x4d 0x54 0x68 0x64. If the file loaded is a MIDI file, the parsing process starts. Otherwise, a message is sent to the Papyrus engine to terminate the parsing process (F2. B).

The MIDI parser begins by loading the chunk header hexadecimal bytes and calling the function that gets the MIDI type, the number of tracks, and the PPQN of the loaded file (F2. C). All those information are sent to the Papyrus modeler, which in turn saves them to an external database. If the number of tracks is higher than one, a message is sent to Papyrus to terminate the process, since the developed system is expecting a one-track melody. After receiving the chunk header data, the system begins reading the next bytes. If the next byte is equal to 0xff, the system recognizes that the following sequence of bytes might be the information associated with the tempo, the time signature, or the key signature. If the next byte is equal to 0x51, the system loads the next four bytes, calculates the tempo, then sends the tempo of the song to Papyrus. If the next byte is 0x58, the system loads the next five bytes, analyzes them, then sends the time signature to Papyrus (F2. E). If the next byte is 0x59, the system gets the key signature from the next three bytes, then sends the key signature used in the MIDI file. After retrieving the essential information, the system starts reading the notes and sending them sequentially to Papyrus to be saved (F2. D). Each note has a timestamp called delta time that changes according to the assigned tempo. A turned-on note (0x90) keeps sounding until receiving a Note OFF message (0x80). Both bytes (0x90 and 0x80) are followed by two bytes that indicate the note and its velocity. At the end of the track (0xff 0x2f 0x00), the parsing system sends a message to the Papyrus engine indicating that the parsing process has been successfully completed (F2. F).

Bars (or Measures) Decomposition Process

In order to split the music into equal measures, from rhythmic perspectives, one should know the numerator and denominator of the time signature for the sequence of consecutive measures. As shown in Figure 4, the duration of a whole note is equivalent to the duration of two half notes, the duration of two half notes is equivalent to the duration of four quarter notes, etc. The time signature is a fraction representing the quantity and the value of the beat. For example, a 4 4 (or 4 × 1 4 ) time signature (Figure 3) indicates that we have 4 beats per bar, each one of them having the value of one quarter note (1/4 of a whole note). A 2 2 (or 2 × 1 2 ) time signature indicates that we have 2 beats per bar, each one of them having the value of one-half note (1/2 of a whole note). To facilitate the development of the bars decomposition system, we have normalized the calculation so that the number of quarter note N per bar is obtained, whatever the time signature is. The decomposition process places the notes into numbered bars. It will not affect the musical composition, the notes, and the arrangement process. It only assembles the notes according to their values and the time signature. Considering that n is the numerator and d is the denominator of the time signature, the number of quarter notes per measure is calculated as follows: 𝑁 = 4 × 𝑛 𝑑 For example, in a the total number of existing notes (TN) (F4. B), then calculates the number of quarter notes (N) required for each bar. Afterwards, the system initiates two variables c1 and c2. While scanning the notes of the melody sequentially and based on the time signature, the variable c1 counts the number of equivalent quarter notes required to fill the bars. The counter c1 is reset to 0 each time it reaches the value N. The variable c2 is in charge of counting the total number of scanned notes. During the scanning process, if the counter c2 is lower than TN and the counter c1 reaches N, a new bar is generated, and then filled in again. This process is repeated until the end of the track (F4. C, F4. D, F4. E).

Three main challenges that should be taken into consideration are the following ones: (1) the first bar of a melody starts in anacrusis (known as "levare" in Italian): in this case, the first bar in anacrusis must be eliminated and the arrangement process starts from the beginning of the second bar. The anacrusis simply entails of the unaccented note/notes preceding the first accent of the rhythmic division in a music composition; (2) the duration of the last bar notes is less than the expected duration: in this case, musical rests must be added according to the remaining unfilled duration; (3) pair tied notes exceeding the remaining rhythmic value available in the bar: in this case, the note is duplicated, then assigned the corresponding duration in the current bar and the next new established bar.

Scale Detection

The scale detection stage is considered as the most delicate phase in the automated music arranger system, since incorrect scale detection leads to incorrect chord generation, and subsequently to a bad arrangement result. A scale detection component is developed to detect the six scale possibilities based on the last played notes in the MIDI file. Knowing the last note (n) of the melody, this component gets the last second, third, fourth, fifth, sixth, and seventh intervals in relation to the note (n), checks if the fourth and fifth are perfect intervals or augmented, then checks the other intervals if they are minor or major. Based on the results of the aforementioned intervals, the system assigns the right scale to the uploaded melody (see Figure 5, that describes the BPMN model of the scale detection component).

Based on the last note of the melody and the scanned intervals, six possibilities exist; i.e., if the last note of the melody is C, the melody can be on: (1) C major if the melody ends with the tonic; (2) C minor if the melody ends with the tonic; (3) A minor if the melody ends with the minor third of the tonic; (4) F major if the melody ends with the right fifth of the tonic; (5) F minor if the melody ends with the right fifth of the tonic; or (6) Ab major: if the melody ends with the major third of the tonic.

The Simulation Process for the Chord Progression Selection

After choosing the scale of the melody, different chord progressions are simulated, analyzed, and then selected. In this section, we define the system component in charge of simulating all possible chord progressions to finally select the best chord progressions that accompany the melody. Some chord progressions are built for major melodies and others for minor melodies. The main advantage of an arrangement based on well-defined chord progressions is that it does not produce a list of successive chords that leads to poor music arrangement results. As the number of chord progressions increases, the repetition of chords decreases. This component assigns to each bar one or more chords. The melody usually starts with an accompaniment of I_Maj or I_min chord, depending on the scale. "I" means the tonic note, "Maj" for major chords, and "min" for minor chords. The system then proceeds to the simulation and comparison of major chord progressions for major scales, and minor progressions for minor scales. Chord progressions are of different sizes. They can be composed of 2 chords up to 64 consecutive chords. Hundreds of possible progressions for minor and major melodies exists. The chord progressions are saved in an external database linked to Papyrus. All progressions are chosen based on thorough theoretical and harmonic studies. The system scans the beats/bars of the melody successively and assigns one or multiple chords to each beat/bar. The simulated first beat/bar of the melody is compared with the first chord of each progression. Next beat/bar is compared with the second chord of each progression, and so on. For each simulated chord, a level of reliability (LRC) is assigned to it. The chord that has a LRC lower than the required level is rejected with the progression that contains it. A chord progression is considered acceptable when all its chords have an LRC higher than the required level. All selected progressions (not previously rejected) are considered eligible to accompany the melody. The same set of bars can have more than one accepted progression. Therefore, another analysis takes place to choose the best chord progression between the selected ones. By calculating the mean value of LRCs, one can have the level of acceptance for each chord progression (LAS). Based on the LAS, the system chooses the most appropriate chord progression that accompanies a part of the melody. When a progression is chosen, the whole process starts again until the end of the track is reached (see BPMN model in Figure 6). Furthermore, different factors and parameters affect the selection of chords/progressions: (1) the length of chord progressions: longer chord progressions (having a greater number of chords) are preferred, i.e., a 6-chord progression having the same LAS as a 2-chord progression will be selected; (2) when common notes are found between the notes of the chord and those of the corresponding beat/bar, the LRC is increased; (3) the LRC is decreased when there are notes in the chord dissonant with the notes of the beat/bar. For example, minor, and sometimes major, second notes; and (4) a chord repetition can occur between consecutive beats. However, the repetition should be avoided between two consecutive bars. Once the chord progressions are chosen successively, and the style of music is determined on Papyrus, the system starts calling the respective MIDI arrangements; e.g., if the user chooses the "Samba" style, the system retrieves all required Samba chord-based arrangements saved in advance, offline. The arranged chords are loaded and concatenated in a sequential order, based on the chord progressions simulation result.

Figure 6: Chord progressions simulation, analysis, and selection process.

Papyrus Modeler Orchestration Process

The BPMN modeling language has been used to build the different components that compose the complex automated music arrangement system. However, the implementation of this system needs an interconnection interface to ensure the communication between several types of processes and technologies presented in the aforementioned sections. For this purpose, we use Papyrus (Guermazi et al. 2015) , an open-source project, providing an integrated environment for editing and executing processes through Unified Modeling Language (UML) and Foundational Subset for Executable UML Models (fUML). The tool offers many interconnection capabilities with complex and heterogeneous systems [START_REF] Gorecki | A Multicomponent Distributed Framework for Smart Production System Modeling and Simulation[END_REF]2021). Papyrus supports UML profiling mechanisms allowing one to extend metaclasses of existing UML metamodels by using stereotypes, so to create a so-called UML profile, and to define a custom behavior at the simulation time using the fUML execution engine (Moka). This feature allows users to define a contextdedicated modeling language, and an adapted execution process. For the orchestration process of the music system components, UML profiles have been defined to store the exchanged data between the Papyrus model, considered as the master component, and each of the developed components responsible for the arrangement process. Papyrus adds new parameters to the executed tasks. Several stereotypes are applied to the UML actions of the main UML metamodel. Thus, the created profiles can be applied on any task to define a communication interface with the components. Once a task is identified, the exchanged data are stored, as parameters, inside the task itself. The execution engine receives those parameters and sends them back to another process. Once UML profiles are defined, it becomes necessary to develop new extensions for Moka engine to achieve the required interaction with the parameters added through the UML profiles. These extensions are simulated each time an extended task is launched to interconnect with the music developed components.

RESULTS AND DISCUSSION

In this project, we took into consideration the six scale possibilities that can be obtained from the last played note in the melody, unlike the existing automated arrangers that consider only two possibilities from the last played note. The existing applications usually check the last note (n) of the melody, consider it as the tonic note of the melody, then scan the third interval. If it is a minor interval, the automated arranger executes the arrangement as if the melody is on (n) minor. If it is a major interval, the arrangement is executed on (n) major. In the scale detection system component, we consider the six scale possibilities from the last scanned note of the melody, i.e., a melody ending with a note C can be on C major, F major, Ab major, C minor, A minor, and F minor scales even though the highest possibilities are the two scales C major and C minor. After the scale detection process, the system chooses the best chord progressions to use, based on a grading system that gets the highest LRC from the chords, then calculates the best sequences having the highest LAS. Each chord progression has a sequential set of chords used to call the respective MIDI arrangements, which are loaded based on the musical genre selected by the user. After a comparison with existing automated music arrangement tools, the developed system shows multiple improvements starting from the scale detection to the arrangement result. Even with a simple melody, the existing tools, most of the time, showed a misdetection of the scale leading to a dissonant arrangement with the supplied melody, chords repetition that may be possibly repeated up to four successive bars, dissonant chords generation, and bad chord progressions. A detailed comparison between the developed system and the existing tools can be found in [START_REF] Possik | Arrangement Musical Automatique[END_REF].

CONCLUSION AND PERSPECTIVES

The proposed automated arranger has been developed on a model-based collaborative platform integrating UML modeling and execution. The definition of Moka extensions and UML profiles using the Papyrus environment allow the interconnection and orchestration of complex processes. The developed automated music arrangement framework shows various improvements compared to some existing automated music arrangement tools in the market. This is an ongoing project, work is in progress to introduce new approaches and features, such as: (1) adding a reinforcement learning feature for a better configuration of chords and chord progressions; (2) adding a new component that converts the audio format to MIDI, which allows the user to play the melody on an instrument, or sing it on a microphone, instead of providing it as a MIDI file;

(3) adding a new component that converts the resulting arrangement from MIDI to audio format, while allowing the user to choose between local synthesizers, external synthesizers, or virtual studio technology instruments known as VSTi; and (4) allowing the user to modify the simulated chords and progressions through a graphical interface.

Figure 2 :

 2 Figure 2: BPMN model of the MIDI parser.

Figure 3 :

 3 Figure 3: Time signature, note duration, and divisions.

  𝑛𝑜𝑡𝑒𝑠 𝑝𝑒𝑟 𝑏𝑎𝑟The total number of bars is equal to the total number of quarter notes divided by the number of quarter notes per bar.The BPMN of Figure4describes the bar decomposition process. This process is launched by the model orchestrator Papyrus.

Figure 4 :

 4 Figure 4: BPMN model of the bars decomposition process. The main role of this component is to load the notes, analyze them, and place them into the music bars based on the provided time signature. This component loads the time signature from Papyrus (F4. A) and

Figure 5 :

 5 Figure 5: BPMN model of the scale detection component.

Information Systems: State of the Art and Future Challenges." Information Systems and E-Business Management vol. 15 no. 2. https://doi.org/10.1007/s10257-016-0317-8.