Mor Ngom
email: mor.ngom@univ-reims.fr

Daniel Alpay
email: alpay@chapman.edu

Mamadou Mboup
email: mamadou.mboup@univ-reims.fr

Scale-Shift and Harmonic analysis approach to the Mellin transform for Discrete-time signals

Keywords: Discrete-time Scale-Shift, Hyperbolic Blaschke group, Mellin transform, Voice-Transform, Discrete Wavelet Transform

We investigate the scale-shift operator for discrete-time signals via the action of the hyperbolic Blaschke group. Practical implementation issues are discussed and given for any arbitrary scale, in the framework of very classical discrete-time linear filtering. Our group theoretical standpoint leads to a purely harmonic analysis definition of the Mellin transform for discrete-time signals. Explicit analytical expressions of the atoms of the discrete-time Fourier-Mellin decomposition are provided along with a simple algorithm for their computation. The so-defined scale-shift operator also allows us to establish a mathematical equivalence in between the discrete-time wavelet coefficients of a given discrete-time signal and the corresponding Voice-transform generated by a well-chosen unitary representation of the Hyperbolic Blaschke group, in the classical Hardy space of the unit disc.

Introduction

Shift operators are very basic ingredients of many fundamental concepts in signal processing, such as periodicity, stationarity, linear filtering, etc. These operators also play a key role in several chapters of mathematics including functional analysis and harmonic analysis which fuel the development of linear systems theory and provide most of the signal processing tools for signal representation, signal analysis and signal synthesis.

The two most common shifts, namely translation (lag operator) and scaling, are combined in the definition of the wavelet transform [START_REF] Meyer | Wavelets and operators[END_REF][START_REF] Daubechies | Ten lectures on wavelets[END_REF]. These two shift operators are linked by an anamorphosis, called the Lamperti transformation [START_REF] Lamperti | Semi-stable stochastic processes[END_REF], which corresponds to a logarithmic deformation of the time axis (see also [START_REF] Cohen | The scale representation[END_REF] where the logarithmic distorsion appears as well in the defined scale operator). This anamorphosis defines a map which allows one to express almost any concept related to time-shift to its scale-shift counterpart. For instance, the Mellin transform appears as the image of the classical Fourier transform and the self-similarity property can be derived from the classical stationarity, up to some multiplicative factor [START_REF] Yazıcı | A class of second-order stationary self-similar processes for 1/f phenomena[END_REF], [START_REF] Mboup | On the structure of self-similar systems: A Hilbert space approach[END_REF], [START_REF] Borgnat | Scale invariances and Lamperti transformations for stochastic processes[END_REF]. Unfortunately, this map applies only in the continuous-time setting. Although scaling is used in many applications, the scale-shift operator does not admit a unanimous definition in the discrete-time case. In this case, scaling is generally restricted to some particular and favourable setting, as for instance the dyadic scales in discrete wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] or approximated, as with the principle of aggregation in the study of self-similarity [START_REF] Leland | On the self-similar nature of Ethernet traffic (extended version)[END_REF]. Yet, a clear-cut mathematical definition was given in [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF] in terms of the action of an Abelian subgroup of the group of automorphisms of the unit disc [START_REF] Ford | [END_REF] (also called the hyperbolic Blaschke group). Below, we use this definition to revisit the wavelet transform in the discrete-time setting.

However, we mention that the same operation has already been used for the purpose of frequency warping. Indeed, scaling in time also corresponds to scaling in frequency. Thus, the wavelet transform in discrete time has been revisited by [START_REF] Evangelista | Discrete frequency-wraped wavelets: theory and applications[END_REF], [START_REF] Evangelista | Audio effects based on biorthogonal timevarying frequency warping[END_REF] in the frequency domain, using the frequency warping in associ-ation with filterbanks. Surprisingly, only the frequency warping effect of the operator has been considered so far while the interpretation of the scale-shift in the time domain did not seem to hold the interest of the authors of [START_REF] Evangelista | Discrete frequency-wraped wavelets: theory and applications[END_REF], [START_REF] Evangelista | Audio effects based on biorthogonal timevarying frequency warping[END_REF] (see also [START_REF] Hsu | Flexible dynamic time warping for time series classification[END_REF] and [START_REF] Schmidt | Iterative twodimensional signal warping-Towards a generalized approach for adaption 28 of one-dimensional signals[END_REF]). Now, the group action viewpoint of the scale-shift enlarges the perspective (see [START_REF] Alpay | A characterization of Schur multipliers between character-automorphic Hardy spaces[END_REF] and [START_REF] Alpay | Discrete-time multi-scale systems[END_REF]) and leads to an exact and purely harmonic analysis approach of the Mellin transform for discrete-time signals.

The proposed approach is in contrast with the principle of discretizing integral operators stemming from the continuous-time setting, using some adequat (geometric) grids (see e.g. [START_REF] Bertrand | Discrete Mellin transform for signal analysis[END_REF], [START_REF] Sena | A fast mellin and scale transform[END_REF]). The presented viewpoint also leads us to consider the unitary representations of the hyperbolic Blaschke group and to find that the so-defined discrete-time scale-shift operator is the very basic ingredient of the Voice transform for this group (see Definition 4.1 and [20], [START_REF] Pap | Properties of the voice transform of the Blaschke group and connections with atomic decomposition results in the weighted Bergman spaces[END_REF]).

Note that strong analogies between the wavelet transform and the Voice transform for the hyperbolic Blaschke group have already been observed by several studies [START_REF] Soumelidis | Applying hyperbolic wavelet constructions in the identification of signals and systems[END_REF], [START_REF] Pap | Hyperbolic wavelets and multiresolution in H 2 (T)[END_REF]. We make this link explicit. Indeed, Theorem 4.2 establishes that the discrete wavelet transform is mathematically equivalent to the Voice transform for an Abelian subgroup of the hyperbolic Blaschke group.

Besides the introduction, this paper is divided into four sections. In section 2, we describe briefly the long path leading to the rigorous definition of the scale-shift operator on discrete-time signals then, we give its implementation and illustrate it with some examples. The Mellin transform for a purely discrete-time signal is defined in section 3. Explicit analytical expressions are given for the atoms of the transform, along with a simple algorithm to compute them recursively. Section 4 starts by basically lightening up the voice transform theory. Then, the equivalence between discrete time wavelets and voice transform is established. We conclude our development in section 5.

The scale-shift in discrete-time

The scale-shift operator for discrete-time signals

In discrete-time, scaling has been first well defined in [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF] using a parallelism with the continuous-time case and an algebraic viewpoint. Without loss of generality, all signals considered in this work are real valued. To begin, recall that, given a signal f (t), t ∈ R, its scale-shift at scale α ∈ R + * reads as f (αt). Also, recall that if

F (s) is the Laplace transform (LT) of f (t), t ∈ R + , then the LT of its scale-shifted version f (αt), α > 0, is 1 α F s α .
Scaling is thus expressed, via the operator S α : u → S α (u) = αu, in the same way in the two domains of continuous-time (u ∈ R + ) and frequency (u ∈ C + ).

In discrete-time world however, scaling a given signal f = {f n } n 0 does not come directly as it is in continuous-time. The key step leading to its clear mathematical definition, given in [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF], was the observation that the operator S α is equivalent to the action of a subgroup of the hyperbolic Blaschke group on the unit disc D of the complex plane:

G = γ(z) = γ 1 z + γ 2 γ2 z + γ1 , |γ 1 | 2 -|γ 2 | 2 = 1 (1)
Each element of G maps D, the unit circle T and the exterior of the closed

disc D = D ∪ T into itself respectively. For f = {f n } n 0 ∈ 2 , consider its Z-transform F (z) = n 0 f n z n , (2)
which converges to an element of the Hardy space of the disc H 2 (D). Now, for a given fixed θ with |θ| < π 2 , it is well known that the Möbius function G θ (s) = e iθ -s e -iθ +s maps C + into D conformally [START_REF] Ford | [END_REF]. Then, for any z ∈ D, there exists a unique s ∈ C + such that z = G θ (s). Hence, from F (z), one considers the map defined in C + by

X(s) := F (G θ (s)) (3)
The map X is well defined and it coincides with the LT of some causal continuoustime signal x(t), t 0 (this is by one of the fundamental results in Fourier-Laplace analysis [START_REF] Stein | Fourier analysis: an introduction[END_REF]). In the frequency domain, the scale-shifted signal x t α , α > 0, is represented by αX

(αs) = X (S α (s)) = F (G θ • S α (s)).
Replacing s by G -1 θ (z) gives

X S α • G -1 θ (z) = F G θ • S α • G -1 θ (z) (4)
This establishes the relation between scaling a discrete-time signal, in the frequency domain, and its continuous-time counterpart. It remains now to identify exactly the discrete-time signal whose frequency representation corresponds to (4) (up to normalization). This leads us to consider the family

G θ = γ α (z) = G θ • S α • G -1 θ (z), |θ| < π 2 . ( 5
)
Now we have:

Lemma 1 ( [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF][START_REF] Alpay | Discrete-time multi-scale systems[END_REF]). Endowed with the composition law •, G θ , forms a group which is a subgroup G. Each γ α ∈ G θ can be written as

γ α (z) = G θ • S α • G -1 θ (z) = γ 1 z + γ 2 γ2 z + γ1 , ( 6
)
where

γ 1 = e iθ + αe -iθ and γ 2 = e iθ (1 -α). ( 7
)
Conversely, every element γ(z) ∈ G is in this form for some θ = θ γ and some

α = α γ .
Fixing θ makes G θ an Abelian group. In all the sequel we fix θ, and we note α γ for the scale associated by [START_REF] Mboup | On the structure of self-similar systems: A Hilbert space approach[END_REF] to the element γ of G θ (also called the multiplier of γ, see [START_REF] Ford | [END_REF]). A direct computation shows that for all γ, ϕ ∈ G θ , we have the relation

α γ α ϕ = α γ•ϕ (8)
which establishes a group homomorphism between G θ and the multiplicative group of R + * . From now on, the coefficients of γ α (z) in [START_REF] Borgnat | Scale invariances and Lamperti transformations for stochastic processes[END_REF] are normalized such

that |γ 1 | 2 -|γ 2 | 2 = 1, as in (1)
:

γ 1 = e iθ + αe -iθ 2 √ α cos θ (9) γ 2 = e iθ (1 -α) 2 √ α cos θ . ( 10
)
Following [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF], we introduce the Hardy space unitary operator T αγ : H 2 (D) → H 2 (D) defined by

T αγ F (z) := F (γ(z)) γ2 z + γ1 = n 0 f n Φ (γ) n (z) (11)
where

Φ (γ) n (z) = [γ(z)] n γ2 z + γ1 , n = 0, 1, . . . ( 12
)
The transform [START_REF] Ford | [END_REF] is just a normalization of (4). For each n 0, it can be readily seen that Φ (γ) n (z) is the image of z n under T αγ . Being unitary [START_REF] Alpay | Discrete-time multi-scale systems[END_REF], T αγ preserves the scalar product and therefore the family {Φ

T αγ F (z) = n 0 f (γ) n z n ( 13
)
where ∀n 0,

f (γ) n = z n , T αγ F (z) H 2 (D) = m f m φ (γ) n,m , (14)
with

φ (γ) n,m = z m , Φ (γ) n (z) H 2 (D) . ( 15
)
Therefore, we have the following Definition. The discrete-time scale shift operator is defined as the mapping

f = {f n } n 0 -→ f (γ) = {f (γ) n } n 0,γ∈G θ . ( 16
)
described in (11)- [START_REF] Schmidt | Iterative twodimensional signal warping-Towards a generalized approach for adaption 28 of one-dimensional signals[END_REF]. In matrix form, this reads as

f (γ) = Φ (γ) f ( 17
)
where the matrix Φ (γ) = φ This operator allows one to make a direct zoom in any given discrete-time signal f = {f n } n 0 for any given scale 1 αγ . Now that the role of G θ in this operator is clear, we invoque the homomorphism ( 8), ( 9) and [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF] to replace all functions of γ ∈ G θ by functions of α ∈ R + * . Thus, if {f n } is the discrete-time analog of f (t) then we will denote by {f

(α) n } the discrete-time analog of f (t/α), keeping in mind that α = α γ for some γ(z) ∈ G θ as in Lemma 1.

Implementation of the discrete-time scale shift

The implementation is based on a digital filtering process. Note that a related implementation appeared in [START_REF] Heuberger | A generalized orthonormal basis for linear dynamical systems[END_REF] in the context of system identification with Laguerre basis (see also [START_REF] Manngård | Identification of low-order models using Laguerre basis function expansions[END_REF] and [START_REF] Wahlberg | On approximation of stable linear dynamical systems using Laguerre and Kautz functions[END_REF]). We start by bringing back the functions Φ (γ) n (z) in [START_REF] Evangelista | Discrete frequency-wraped wavelets: theory and applications[END_REF] and observe that they satisfy the recurrence relation

Φ (γ) 0 (z) = 1 γ2 z + γ1 (18)
Φ (γ) n (z) = γ(z)Φ (γ) n-1 (z), n 1. ( 19
)
Let {δ n } denote the unit impulse sequence. Given a discrete-time signal f = {f n }, its scale-shift at a scale 1 αγ > 0, as given in [START_REF] Hsu | Flexible dynamic time warping for time series classification[END_REF], can be implemented by the scheme below: 11), ( 18) and ( 19)

δ n Φ 0,γ (z) φ n,0 γ(z) φ n,1 γ(z) • • • φ n,2 • • • ⊗ f 0 ⊗ f 1 ⊗ f 2 ⊕ ⊕ • • • • • • ⊕→ f (αγ ) n Figure 1: Implantation diagram obtained from (

Finite support and Gibbs phenomenon

Without loss of generality, we consider the example below to fix the ideas.

Let {f n } be the discrete-time signal defined for all n 0 by f n = f (nT s ) where T s = 0.01 and

f (t) = 0.66 (t -1.87) 2 + 0.1 e -0.66(t-2) 2 . ( 20
)
Below, we consider the signal {f n } N n=0 and investigate the effect of the finite support N on its scale-shitfs f

(α) n Nα n=0
, for various N and α. As one may readily observe, the support of the scale shifted signal is given by N α = N α , where the notation • stands for the integer part of the argument.

Smooth signal

In this first example, we set N = 400. This is large enough to consider that the whole support of f (t) is covered. The top plot in figure 2 } N n=0 defined by

h n = f (α) n (21) e (α)•(β) n = f n -h (β) n ( 22
) with α = 1 β = 1 √ 2
for the red curve and α = 1 β = √ 2 for the blue curve. The error signal e

( 1 √ 2 )•( √ 2) n n 0
is negligible, which thus shows an almost exact reversibility of the discrete-time scale-shift operation. The same conclusion can be drawn for the other error signal, e

( √ 2)•( 1 √ 2 ) n n 0
, although it shows some fluctuations with gradual amplitudes in the second half of its support. As explained below, these fluctuations come from the Gibbs phenomenon, induced by the fact that f N > 0 corresponds to a jump, even though the amplitude is small.

Nonsmooth signals and Gibbs oscillations

Recall that, as a transfer function, γ(z) in figure 1 has an infinite impulse response. Strictly speeking, the matrix Φ (γ) in ( 17) is therefore doubly infinite even though the implementation considers a limited row dimension of N α as given above, to keep only the meaningful part of the scaled signal. The α-scale shift of a finite support signal {f n } N n=0 then reads in matrix form as

               f (α) 0 . . . f (α) Nα * * . . .                = Φ (γ)                   f 0 . . . f N 0 0 0 . . .                   (23)
where the " * " are meaningless values since only the first N α + 1 components are of interest. This equation shows that the scale shift operator sees a jump in the finite support original signal, by considering its extension of infinite duration by zero padding. The jump is always present unless the last samples smoothly decay to zero, that is: f n → 0, n = N 0 , N 0 + 1, . . . , N for some range N 0 < N . As is well known, low-pass filtering a signal with jump induces a Gibbs phenomenon. Now, we claim that these Gibbs oscillations start to occur, in the scaled version of the signal, from the time index

n α = N α . ( 24
)
When α < 1 (scale contraction) we have n α > N α and thus the Gibbs phenomenon does not appear in the finite support of the scaled signal. However, it happens in the meaningless part of the scaled signal, marked by the symbol " * " in the left-hand side of [START_REF] Pap | Hyperbolic wavelets and multiresolution in H 2 (T)[END_REF].

For α > 1 however, the oscillations become visible in the whole time index

interval [n α , N α ].
These conclusions are illustrated in the next example.

In this second example, we reconsider the discrete-time signal {f n } N n=0 sampled from f (t) in [START_REF] Pap | The voice transform on the Blaschke group II[END_REF], with a finite duration N such that f ((N + m)T s ) is significantly different from 0, for all m = 0, 1, . . . , M , for some M > 0. For example, N = 345 satisfies this condition, with M = 55. The jump in the original signal, displayed by the black colored curve in figure 3, induces Gibbs oscillations for all scales α 1 = √ 3 and α 2 = 1 √ 3 . With α 1 > 1, the oscillations occur in the interval [n α1 , N α1 ], as claimed above (see the blue colored curve in figure 3). In the red colored curve of the same figure, corresponding to α 2 < 1, the Gibbs phenomenon is still present but the oscillations start after the time index N α2 since in this experiment, we have n α2 > N α2 . Observe also that as α 1 = 1/α 2 , we have a coincidence of n α1 and N α2 . These results confirm the previous analysis on the effect of a jump on the scale operator. Note also that the analysis already suggests a solution to make the scale shift operator safe from the Gibbs phenomenon.

Mitigating the Gibbs phenomenon

To begin, recall that for α < 1, the Gibbs oscillations are unseen because they happen after the time N α , that is in the meaningless part of the scaled signal (see equation ( 23)). Now, observe that one can get back to the same situation (with jump at the end, as explained before), with α > 1. First we define the

signal { f n } N n=0 of duration N = N + K α such that f n =      f n , n = 0, 1, • • • , N f N n = N + 1, • • • N = N + K α ( 25
)
where K α is some constant depending on α. According to eq.( 23) and setting

N α = N α , the α-scale shift of { f n } N n=0 then reads as                         f (α) 0 . . . f (α) Nα f (α) Nα+1 . . . f (α) Nα * . . .                         = Φ (γ)                       f 0 . . . f N f N . . . f N 0 . . .                                  K α lines (26)
We thus get back to the setting in eq.( 23) where the original jump, at time index N , is now carried over to the time index N = N + K α . Accordingly, the Gibbs oscillations start to occur at the time index n α = N α .

N α = N α n α = N + K α α .
Now this is achieved when

K α N (α 2 -1) .
This is confirmed by the results in figure 5, where the scaled versions {f

(α) n } Nα n=0
of the finite support signal {f n } N n=0 (black curve) are free from any Gibbs oscillations for all α, be it lower or greater than 1. The scaled version is obtained from the filtering scheme in figure 1 equation ( 23), for α < 1 (red curve in the plot). For α > 1, it is defined by f

(α) n = f (α) n , n = 0, • • • , N α (blue curve).

Mellin transform for discrete-time signals

To outline the guiding ideas of our approach, we begin with a brief review of the basic elements of harmonic analysis that underlie the abstract Fourier transform. This subject is standard and is exposed in a number of classical (re-edited) textbooks such as e.g. [START_REF] Rudin | Fourier analysis on groups[END_REF], [START_REF] Goldberg | Fourier transforms[END_REF] and [START_REF] Hewitt | Abstract Harmonic Analysis[END_REF]. Nonetheless, we recall the basic bricks because the subject may not be familiar to some readers. Next, we

show how easy and natural the definition of the Mellin transform for discretetime signals falls within this general framework. The last part of this section will provide a simple algorithm to compute the atoms of the decomposition.

Abstract Harmonic Analysis

Basic definitions

Let (G , •) be a locally compact Abelian group, with Haar measure µ and

neutral element ι. A function σ : G → T satisfying i) σ(ι) = 1 (27) ii) σ(a • b) = σ(a)σ(b), ∀a, b ∈ G (28)
is called a character. The set of all continuous characters forms a group G called the dual group of G . By the Pontryagin duality, G is also a locally compact Abelian group and we denote by μ its associated Haar measure. Now, for f : G → C in L 1 (G , dµ), the Fourier transform F and its inverse F -1 are defined respectively by

f (σ) = F{f }(σ) = G f (a)σ(a)dµ(a), σ ∈ G (29) f (a) = F -1 { f }(a) = G f (σ)σ(a)dμ(σ), a ∈ G . ( 30
)

Group action

Consider the action of G on a non-empty set E. This can be described via for each different t. In some situations however, it becomes straightforward to determine h t (σ) based on the knowledge of h τ (σ) for any given pair (τ, t) ∈ E 2 .

a morphism η : G × E → E such that • η(ι, t) = t, for all t ∈ E • η(a, η(b, t)) = η(a • b, t),
This arises when the group action is transitive. Recall that the action of G on E is said transitive if for any pair (τ, t) ∈ E 2 , there exists a ∈ G such that τ = η(a, t). Let then t 0 ∈ E be fixed once for all and given t ∈ E, let b t ∈ G be defined such that η(b t , t 0 ) = t. Then we get

h t (a) = h(a, t) = f (η(a, t)) = f (η(a, η(b t , t 0 ))) = f (η(a • b t , t 0 )) = h t0 (a • b t ),
showing that the two variables merge into one and the function h(•, •) reduces to a function of a single variable. The Fourier transform of h t (•) becomes:

h t (σ) = G h t (a)σ(a)dµ(a) = G h t0 (a • b t )σ(a)dµ(a) = G h t0 (a)σ(a • b -1 t )dµ(a • b -1 t ) = σ(b t ) h t0 (σ),
where the last equality is obtained using ( 27)-( 28) and the translation-invariance property of the Haar measure. An example of such simple situation is the case when G is the additive group of real (resp. integer) numbers, and E = G :

The Haar measure is the Lebesgue measure (resp. counting measure) and the characters are the exponential functions σ ω (a) = e -iωa . This corresponds to the classical Fourier transform. Another common example is obtained with the multiplicative group of positive real numbers, with again, E = G : The Haar measure is dµ(a) = da a , the characters are the functions σ ω (a) = a -iω , ω ∈ R and we recover the Mellin transform.

Scaling and Mellin transform in discrete-time

As observed already, scale translation in discrete-time is expressed via the action of the group G θ (θ being fixed as before). We recall that by Lemma 1,[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and Definition 2.1, this action is equivalent to that of the multiplicative group of positive reals, with E = N. The associated morphism η(α, n) is defined such that for a given discrete-time signal f = {f n } n 0 with its Z-transform

F (z) = f n z n , we have f η(α,n) := f (α) n = z n , T α F (z) H 2 (D) . ( 31
)
Clearly, this action is not transitive. Therefore, the two variables α ∈ R + * and n ∈ N cannot be merged.

f m (ω) = M m {f }(ω) = n 0 f n E m,n (ω), ( 32
)
where the n th atom of the transform is

E m,n (ω) = (-1) m 2π cosh(πω) P m (iω)P n (iω), n = 0, 1, 2, • • • ( 33
)
where P n (s) is the polynomial of degree n which reads as

P n (s) = 1 n! n k=0 n k n =1 ( -k - 1 2 -s) (34)
2) -The transform is reversible and the original signal {f n } can be retrieved by

f n = Φ (α) n (z), F (α) (z) H 2 (D) , ( 35
)
where

F (α) (z) = m M -1 f m (s) z m
and where we use the notation M -1 for the classical inverse Mellin transform.

Proof. The proof is given in the appendix.

In the sequel, we define the polynomials

p k,n (s) = n =1 ( -k -s - 1 2 ), n ∈ N, 0 k n (36)
so that the polynomial in (34) reads as

P n (s) = 1 n! n k=0 n k p k,n (s).

Computation algorithm

The following lemma shows some relations on these polynomials that are useful for the computation of the atoms E m,n (s).

Lemma 2. The polynomials p k,n satisfy the following recurrence relations

p n-k,n (s) = (-1) n p k,n (-s) (37) p k,n (s) = 2(n -k) -1 2 -s p k,n-1 (s), ∀n 1 (38) p k,n (s) = -2k+1 2 -s 2(n-k)+1 2 -s p k-1,n (s), 1 k n (39)
Proof. See the appendix.

Lemma 3. On the imaginary axis, s = iω, ω ∈ R, the polynomial P n (s) is either real valued or purely imaginary and we have:

1. For n even, n = 2q,

P 2q (iω) = 2 (2q)! q-1 k=0 2q k {p k,2q (iω)} + 1 (q!) 2 p q,2q (iω) (40)
and p q,2q (iω) is real.

2.

For n odd, n = 2q + 1,

P 2q+1 (iω) = 2i (2q + 1)! q k=0 2q + 1 k {p k,2q+1 (iω)} (41)
These lemmas show how one can compute the polynomials P n (s) for all n, using simple recurrences iniatialized from p 0,0 (s) = 1. Indeed, for each

given n > 0, the polynomials p k,n (s), k = 1, • • • , n can be easily obtained recursively from p 0,n (s) by following the recurrence relation (39). Now, the polynomial p 0,n (s), in turn, easily derives from p 0,0 (s) by a recursive application of the relation (38). Finally, the polynomials p k,n (s), and hence P n (s), can be computed using the diagram in figure 6 below, where the vertical arrow with the circle mark in the middle represents the relation (39) and the one with the square mark represents the relation (38). Based on Lemma 3 and (33), it holds that E m,n (ω) is real when n+m is even and purely imaginary otherwise. The supports and the shapes of the different atoms vary markedly with m and n. This is illustrated in figure 7 below, where the atoms or their imaginary part where appropriate, are represented for some values of n and m.

p 0,0 | P 0 . . . p 0,2(q-1) • • • p q-1,2(q-1)| P 2(q-1) p 0,2q-1 • • • p q-1,2q-1 | P 2q-1 p 0,2q • • • p q-1,2q p q,2q | P 2q p 0,2q+1 • • • p q-1,2q+1 p q,2q+1 | P 2q+1 . . . . . . . . .

Voice transform and wavelet coefficients

In this section, we use the scale-shift operator defined above to revisit two popular signal transforms: the discrete wavelet transform and the Voice transform for the Blaschke group. We establish a strict equivalence between these two transform.

The Voice transform

The Voice transform was originally introduced in [START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decompositions, I[END_REF][START_REF] Feichtinger | Banach spaces related to integrable group representations and their atomic decompositions[END_REF] and serves as a unifying theory for lot of atomic representations (see also [START_REF] Pap | Properties of the voice transform of the Blaschke group and connections with atomic decomposition results in the weighted Bergman spaces[END_REF][START_REF] Pap | The voice transform on the Blaschke group II[END_REF][START_REF] Soumelidis | Applying hyperbolic wavelet constructions in the identification of signals and systems[END_REF]). Generally speaking, one starts by representing a given locally compact topological group G in a given Hilbert space H by unitary operators U g , g ∈ G , defined in H. More precisely, the voice transform of an element f ∈ H generated by U g , g ∈ G , and with respect to some parameter h ∈ H, is defined as

[V h f ](g) = f, U g h H ( 42
)
where •, • H is the inner product in the Hilbert space H.

Related to the present context, let us consider the subgroup G θ defined in [START_REF] Yazıcı | A class of second-order stationary self-similar processes for 1/f phenomena[END_REF], with an angle θ ∈ (-π/2, π/2) fixed as before. Recall that the operators [START_REF] Ford | [END_REF] are unitary representations of G θ in the Hilbert space of the disc, H 2 (D).

T αγ : H 2 (D) → H 2 (D) defined in
Therefore, we can define the Voice transform with respect to G θ , for discretetime signals:

Definition. Let X(z) = n 0 x n z n ∈ H 2 (
D) be given. For a causal discretetime signal f = {f n } n 0 ∈ 2 , we define its associated Voice transform, generated by T αγ , with parameter X, as:

[V X(z) f ] : G θ → C, such that [V X(z) f ](γ) = F (z), T αγ X(z) H 2 (D) = n f n x (αγ ) n ( 43
)
where F (z) = n 0 f n z n

Voice transform to Wavelets coefficients

We are now ready to state and prove the mathematical equivalence between the Voice transform and the wavelet coefficients of a given discrete-time signal.

We mention that the relation in between these two transforms has been studied in [START_REF] Soumelidis | Applying hyperbolic wavelet constructions in the identification of signals and systems[END_REF] without touching an explicit and strict equivalence. The reason lies probably on the lack of a proper scale-shift operator. Indeed, because of the non commutativity of the whole Blaschke group G, different elements of the group present different rotation angles and this precludes to interpret the action of the entire group in terms of scale-shift. Here, we get around this difficulty by fixing the angle θ and just consider the subgroup G θ of SU (1) which, in turn, is Abelian.

Recall that the wavelet transform of a given continuous-time signal f (t) ∈ L 2 (R) is defined as

[W ψ f ](α, u) = 1 √ α R f (t)ψ t -u α dt, ( 44
)
where ψ is the corresponding mother wavelet, α and u are the scale and timeshift parameters, respectively.

In denote the discrete-time version of ψ t α as defined by the discrete-time scale shift operator in Definition 2.1. Then, the discrete wavelet transform may be readily given by the discrete-time version of eq.( 44) by

c k (α) = n f n ψ (α) n-k . ( 45
)
This expression leads to the following result:

Theorem. The voice transform of a discrete-time signal {f n } n 0 is equivalent to its discrete wavelet transform as given in (45), that is:

c k (α) = V [z k Ψ(z)] f (γ α ) (46)
Proof. Let F (z) = n 0 f n z n and Ψ(z) = n ψ n z n be the Z-transforms of the sequences {f n } n 0 and {ψ n } n 0 respectively. The Z-transform of the discrete-time scale shifted signal ψ (α) n then readily reads as

Ψ (α) (z) = n ψ (α) n z n = T α Ψ(z).
Likewise, for any fixed k ∈ N, we will use the notation Ψ

(α) k (z) = n ψ (α) n-k z n = T α [z k Ψ(z)]
, to represent the discrete-time scale shift of {ψ n-k } n 0 . The proof 20 follows upon noting that the expression (45) also reads in a form of a scalar product

c k (α) = n f n ψ (α) n-k (47) = F (z), Ψ (α) k (z) H 2 (D) (48) = F (z), T α z k Ψ(z) H 2 (D) (49) = V [z k Ψ(z)] f (γ α ), (50)
where γ α (z) is the element of G θ associated with the scale α as in [START_REF] Yazıcı | A class of second-order stationary self-similar processes for 1/f phenomena[END_REF]. This ends the proof.

Conclusion

This The atoms of the decomposition have been given by explicit expressions and a easy computational algorithm has also been devised. Finally, the definition of this operator has enabled us to establish an explicit equivalence between two popular signal transforms: the Voice transform associated to a specific Blaschke hyperbolic subgroup and the discrete wavelet transform.

The corresponding Mellin transform is well-known as

Y z (s) = λ -s z Γ(s)Γ(n + 1 -s) Γ(n + 1) ,
where Γ(•) is the classical Gamma function. Therefore, we deduce the Mellin transform

Φ n (s, z) := M Φ (α) n (z) = 2 1 + z n k=0 (-1) k n k λ k z Y z s + k + 1 2 = 2 1 + z n k=0 (-1) k n k λ k z λ -s-k-1/2 z • Γ(s + k + 1/2)Γ(n -k + 1/2 -s) Γ(n + 1) . Setting s = s + k + 1/2, we write Γ(n -k + 1/2 -s) = Γ(n + 1 -s) = n =1 ( -s)Γ(1 -s)
With the Euler's reflection formula

Γ(ζ)Γ(1 -ζ) = π sin πζ , ∀ζ ∈ C \ Z,
and recalling the expression (36), namely

p k,n (s) = n =1 ( -s k ) = n =1 ( -k -1/2 -s),
we obtain:

Φ n (s, z) = 2λ -s-1/2 z 1 + z n k=0 (-1) k n k π p k,n (s) n! sin [π(s -k -1/2)] = 2λ -s-1/2 z 1 + z π cos (πs) n k=0 n k p k,n (s) n! = 2λ -s-1/2 z 1 + z π n! cos (πs) P n (s) 5.1.2. Computation of z m , Φ n (s, z)
To finish the proof, we need to express the development of Φ n (s, z) in a formal power series of the variable z, s being a fixed parameter. To begin, we set

A s (z) := λ -s-1/2 z 1 + z = (1 -z) -s-1/2 (1 + z) -s+1/2 .
The following formal power series

(1 -z) -s-1/2 = ∞ =0 Γ(s + 1/2) Γ(s -+ 1/2)Γ( + 1) z (1 + z) -s+1/2 = ∞ =0
(-1) Γ(-s + 1/2) Γ(-s -+ 1/2)Γ( + 1) z , allow one to rewrite A s (z) as

A s (z) = ∞ k, =0 (-1) k Γ(-s + 1 2 )Γ(s + 1 2 )z k+ k! ! Γ(-s + 1 2 -k)Γ(s + 1 2 -) = ∞ m=0 m k=0 (-1) k k!(m -k)! Γ(-s + 1 2 ) Γ(-s + 1 2 -k) • Γ(s + 1 2 ) Γ(s + 1 2 -m + k) z m
(-m + k -s - 1 2 ) = n =1 (-+ k -s + 1 2 ) = (-1) n n =1 ( -k + s - 1 2 )
= (-1) n p k,n (-s)

Recurrence on n (38)

p k,n (s) = n =1 ( -k -s - 1 2 ) = n-1 =1 ( -k -s - 1 2 ) • (n -k -s - 1 2 )
= a k,n (s)p k,n-1 (s)

where a k,n (s) = n -k -s -1 2 = 2(n-k)-1 2 -s.

Recurrence on k (39)

p k,n (s) = n =1 ( -k -s - 1 2 )
=

Proof of Lemma 3

Observe first that the symmetry relation (37) implies that p q,2q (iω) = p q,2q (-iω) = p q,2q (iω)

showing that p q,2q (iω) is real. Next using the recurrence relation (37) and the identity n k = n n-k , we can simplify the computation of the polynomials P n (s), for n even/odd, by n! (p k,2q (s) + p k,2q (-s)) + 1 (q!) 2 p q,2q (s). The Lemma follows upon noting that for s = iω, ω ∈ R, it holds that p k,n (-iω) = p k,n (iω) since all coefficients of p k,n are real.

n

  (z)} n 0 is orthogonal and complete in the Hardy space H 2 (D) of the unit disc. Hence, expressing each Φ (γ) n (z) in the standard basis of H 2 (D), we can rewrite T αγ F (z) as

  vectors f and f (γ) contain the samples f n , n 0 and their scale shifted counterparts f (γ) n , respectively.

Figure 2 :

 2 Figure 2: Effectiveness of the scale shift operator even for irrational scales (top plot) -Reversibilty of the scaling (bottom plot)

Figure 3 :

 3 Figure 3: Gibbs oscillations and scale shift of a nonsmooth signal.

with α > 1

 1 by pushing (artificially) the jump of the signal far enough. More precisely, consider the computation of the α-scale shift of the signal {f n } N n=0

Figure 4 : 3 .

 43 Figure 4: Making scale shitf insensitive to Gibbs oscillations

Figure 5 :

 5 Figure 5: Scale-shift immune to Gibbs oscillations

  For any fixed pair (n, m) of positive integers with m = n, f (α) n and f (α) m are therefore two unrelated functions of α. Hence a different Mellin transform needs to be computed for each different m. Theorem. Let f = {f n } n 0 ∈ 2 be a causal discrete-time signal and let m ∈ N be a given time index. Then 1) -The Mellin transform of {f n } n 0 with respect to the time index m, denoted by f m = M m {f }, is given by

Figure 6 :

 6 Figure 6: Computation diagram of the set of polynomials Pn(iω).

Figure 7 :

 7 Figure 7: Some atoms Em,n of the Mellin transform

  paper has investigated a clear cut definition of the scale-shift operator in discrete-time. Its implementation has been given under a very simple and classical linear filtering setup. It has then been illustrated with some examples showing the presence of Gibbs phenomenon under the scale-shift of discretetime signal with finite support. However, the Gibbs flusctuations can be fully avoided and a simple and effective solution is provided. Then, the scale-shift operator has been used to define rigorously the Mellin transform for discrete-time signals. The definition falls within the classical harmonic analysis framework.

=1((

  k m k r k (s)r m-k (-s).with r 0 (s) k by m -k and s by -s, leads to r m-k (-s) = the product in the expression of p k,n (s) in (36) asp k,n (s) = k -k -1/2 -s) n =k+1 -k -1/2 -s)then one immediately recognizes the equalityp k,n (s) = (-1) n-k r k (s)r n-k (-s).This leads to the conclusion thata m (s) = 1 m! m k=0 (-1) k m k (-1) m-k p k,m (s) = (-1) m P m (s).The proof is now complete withE m,n (s) = z m , Φ n (s, z) (54) = π n! cos (πs) P n (s) z m , A s (z)(55)= π n! cos (πs) (-1) m P m (s)P n (s)

- 1 -- 1 ) -s - 1 2 =

 112 (k -1) -sb k,n (s)p k-1,n (s)where b k,n (s) =

For

For n = 2q + 1 ,

 1 k,n (s) + p n-k,n (s)) ,n (s) -p k,n (-s)) .

  the context of this present work, signals are assumed to be known only at discrete-time instants i.e. f = {f n } n 0 . Let {ψ n } be the discrete-time version of the mother wavelet, that may be obtained, for instance, by ψ n = ψ(nT s ) where T s is a sampling time. Accordingly, let ψ

	(α)
	n

The first task towards the computation of M Φ

n (z) as a function of α, regarding n and z as fixed parameters. Without loss of generality, we fix θ = 0 and consider the subgroup G 0 . The dependence of Φ (α) n (z) on α can be made explicit by recalling ( 12), ( 9) and [START_REF] Mboup | A character-automorphic Hardy spaces approach to discretetime scale-invariant systems[END_REF] with θ = 0. We get

where we have set

and where we have used the well-known binomial formula in the last equality.

Next, recall that if Y (s) represents the Mellin transform of a function y(x)

then the Mellin transform of x ν y(x) would read as Y (s + ν). So, let us set


Florence Noble

Dr Nicolas Marie
email: nicolas.marie@parisdescartes.fr

Biased opioid ligands: revolution or evolution ?

Keywords: opioids, biased ligands, analgesia, side effects, TRV130, PZM21

Opioid are the most powerful analgesics ever but their use is still limited by deleterious side effects such as tolerance, dependence and respiratory depression that could eventually lead to a fatal overdose. The opioid crisis, mainly occurring in north America, stimulates research on finding new opioid ligands with reduced side effects. Among them, biased ligands are likely the most promising compounds. We will review some of the latest discovered biased opioid ligands and see if they were able to fulfil these expectations.

Introduction

Opium consumption and its medicinal application dates back to the Neolithic, Bronze, and Iron Age. Ancient Egyptian papyrus records reported the use of opium for pain relief, demonstrating that opioids are used from thousands of years for the treatment of pain. Opium is a substance extracted from Papaver somniferum, and opium poppies produce alkaloids such as morphine [START_REF] Trescot | Opioid pharmacology[END_REF]. In 1973, the notion of "morphine receptors" emerged [START_REF] Pert | Opiate receptor: demonstration in nervous tissue[END_REF], and since that time numerous studies have been performed on these receptors. In the 1990s, opioid receptors were cloned. They are 4 of them, MOR, DOR, KOR and NOR, and are largely found within the central nervous system as well as throughout the peripheral tissues. They are naturally stimulated by endogenous peptides (endorphins, enkephalins, dynorphins and nociceptin/orphanin FQ). All pharmacological studies confirmed that opioid system plays a key role in the modulation of pain perception, but is also involved in many physiological responses, such as respiration, gastrointestinal motility, endocrine and immune functions, which may lead in clinical practice to adverse effects when pain must be managed with opioids.

The study of responses of opioid receptors-deficient mice to opioid agonists clarifies the biological activity of each receptor. The genetic approach unambiguously demonstrated that MOR are essential for all the biological activities elicited by morphine, including for instance analgesia, reward, respiratory depression, constipation, immunosuppression [START_REF] Kieffer | Exploring the opioid system by gene knockout[END_REF]. Thus, a challenge over the past decades has been to try to find MOR opioid ligands able to induce a potent analgesic response, with fewer adverse effects.

MOR belong to the G protein-coupled receptors (GPCR), and it is now well established that once the receptor is activated, several intracellular signaling pathways may be activated. Among these pathways, the focus has been made on two of them, G protein (and its subsequent effectors such as adenylate cyclase and ion channels) and beta-arrestin, and it has been suggested that analgesic effects of MOR ligands are mediated by the downstream Gprotein-dependent pathway, while beta-arrestin-dependent pathway mediates most of the side effects. Thus, it was speculated that development of biased MOR ligands for G-proteindependent pathway may have high analgesic potency, but with fewer undesirable effects. The main goal of this review is to briefly describe this strategy and to discuss its advantages and limits.

Side effects of opioids and interest in biased agonism

With great therapeutic potential, usually comes great side effects. The opiates do not seem to deviate from this rule. Opioid receptors and especially MOR are highly expressed in the brain structures involved in control of breathing including the pre-Bötzinger complex or the Kolliker-Fuse neurons [START_REF] Montandon | PreBötzinger Complex Neurokinin-1 Receptor-Expressing Neurons Mediate Opioid-Induced Respiratory Depression[END_REF]. Activation of MOR in the pre-Bötzinger complex with fentanyl promotes a respiratory depression [START_REF] Montandon | Role of the pre-Bötzinger Complex in opioidinduced respiratory depression in adult rats in-vivo[END_REF] and the knockout of MOR in this complex reduced the effect of morphine on breathing [START_REF] Bachmutsky | Opioids depress breathing through two small brainstem sites[END_REF]. MOR are also expressed in the enteric nervous system and their activation will cause reduction of gastrointestinal motility and thus constipation [START_REF] Wood | Function of opioids in the enteric nervous system[END_REF].

Following a protracted exposure, but also after an acute exposure at high dose, a tolerance to opioid-induced pharmacological effects appears. It noteworthy that the level of tolerance depends on the pharmacological effect studied. For instance, there is evidence to suggest that the tolerance to respiratory depression is weak as compared to the antinociceptive tolerance [START_REF] Mohammed | Comparison of tolerance to morphine-induced respiratory and analgesic effects in mice[END_REF] which might explained the fatal overdose that could occur in opioid users [START_REF] White | Mechanisms of fatal opioid overdose[END_REF]. Finally, a repeated use of opioids, even in the frame of a pain treatment, could lead to development of addictive behaviours. This is well illustrated by the opioid crisis in US where postoperative opioid prescription is a significant contributor to opioid epidemic, and more than 350,000 people died from opioid overdose between 1999 and 2016 (https://www.cdc.gov/drugoverdose/epidemic/index.html). All these side effects are a real impairment in the opioid use, particularly in chronic pain. Among the different strategies to find potent analgesics with limited side effects, the one based on the particular properties of the opioid receptors and their ligands, the biased signalling, has received great attention since few years even though the concept (as known as functional selectivity) is older [START_REF] Kenakin | Biased Receptor Signaling in Drug Discovery[END_REF]. Indeed, the notion of functional selectivity was probably first suggested in a review on serotoninergic receptors where the ability of a single GPCR to active one or more distinct signalling cascades was described [START_REF] Roth | Multiple mechanisms of serotonergic signal transduction[END_REF]. And the term of biased agonist was first introduced by Jarpe and co-workers when they observed that a peptide promoted a biased activity on chemokine receptors [START_REF] Jarpe | d-Arg1,d-Phe5,d-Trp7,9,Leu11]Substance P Acts as a Biased Agonist toward Neuropeptide and Chemokine Receptors[END_REF]. For opioid receptors, the studies of Bohn and co-workers along with this notion opened the way to find new opioid ligands. As rapidly mentioned earlier, following ligand binding, MOR activation could result in the activation of multiple downstream pathways through either G protein dependent processes (e.g. adenylate cyclase inhibition or regulation of ion channels) or G protein independent processes (e.g. beta-arrestin signaling). Beta-arrestins are proteins, existing in two isoforms (beta-arrestin-1 and -2), that bind the activated and phosphorylated receptor and are responsible for its desensitization and internalization [START_REF] Allouche | Opioid receptor desensitization: mechanisms and its link to tolerance[END_REF]. In beta-arrestin-2 knockout mice, Bohn and collaborators found that morphine analgesia was enhanced, with abolition of antinociceptive tolerance [START_REF] Bohn | Enhanced morphine analgesia in mice lacking beta-arrestin 2[END_REF][START_REF] Bohn | Mu-opioid receptor desensitization by beta-arrestin-2 determines morphine tolerance but not dependence[END_REF] and reduction of respiratory depression and acute constipation [START_REF] Raehal | Morphine side effects in beta-arrestin 2 knockout mice[END_REF]. These data suggested that whereas G protein-dependent pathway mediated analgesia, beta-arrestin-dependent signalling pathway promoted side effects. Using different strategies (high-throughput screening, structure-based virtual screening or synthesis-driven approach [START_REF] Bermudez | Strategies for the discovery of biased GPCR ligands[END_REF]), new biased MOR ligands were discovered. They are mostly biased towards G proteins.

Biased ligands are promising compounds

The most notable compound,TRV-130 (Oliceridine) was discovered after a high-throughput screening by Trevena Inc. [START_REF] Chen | Structure-Activity Relationships and Discovery of a G Protein Biased μ Opioid Receptor Ligand[END_REF]. It has a moderate biased toward cAMP pathway (about 3 fold) and was unable to promote MOR endocytosis according to its very weak ability to recruit beta-arrestin-2. Preclinical studies showed that TRV130 had a higher potency than morphine to induce antinociceptive responses with a reduced respiratory depression and gastrointestinal effects, suggesting a better therapeutic index as compared to morphine [START_REF] Dewire | A G protein-biased ligand at the μ-opioid receptor is potently analgesic with reduced gastrointestinal and respiratory dysfunction compared with morphine[END_REF]. Moreover, in the tail immersion test, no tolerance was measured after a three days treatment [START_REF] Altarifi | Effects of acute and repeated treatment with the biased mu opioid receptor agonist TRV130 (oliceridine) on measures of antinociception, gastrointestinal function, and abuse liability in rodents[END_REF]. Clinical studies validated these results. Indeed, in a phase 2 study comparing TRV130 with morphine in patients after abdominoplasty, the biased agonist provided not only a quicker analgesia, but was also associated with fewer side effects including hypoventilation, nausea and vomiting [START_REF] Singla | A randomized, Phase IIb study investigating oliceridine (TRV130), a novel µreceptor G-protein pathway selective (µ-GPS) modulator, for the management of moderate to severe acute pain following abdominoplasty[END_REF]. In a phase 3 study, TRV130 was shown to be efficient in reducing pain intensity in two surgical procedures: bunionectomy (APOLLO-1 trial) [START_REF] Viscusi | APOLLO-1: a randomized placebo and active-controlled phase III study investigating oliceridine (TRV130), a G protein-biased ligand at the µ-opioid receptor, for management of moderate-to-severe acute pain following bunionectomy<[END_REF] and abdominoplasty (APOLLO-2 trial) [START_REF] Singla | APOLLO-2: A Randomized, Placebo and Active-Controlled Phase III Study Investigating Oliceridine (TRV130), a G Protein-Biased Ligand at the μ-Opioid Receptor, for Management of Moderate to Severe Acute Pain Following Abdominoplasty[END_REF]. Analysis of side effects in the pooled cohorts of APOLLO trials showed that whereas TRV130 induced a dose dependent respiratory depression, this effect is lower than the one induced by morphine [START_REF] Ayad | Evaluating the Incidence of Opioid-Induced Respiratory Depression Associated with Oliceridine and Morphine as Measured by the Frequency and Average Cumulative Duration of Dosing Interruption in Patients Treated for Acute Postoperative Pain[END_REF], and that TRV130 had a higher probability of producing analgesia rather than respiratory depression [START_REF] Dahan | Benefit and Risk Evaluation of Biased μ-Receptor Agonist Oliceridine versus Morphine[END_REF]. These analyses also demonstrated a lower risk of experiencing nausea and vomiting in patients treated with TRV130 compared to morphine-treated patients [START_REF] Beard | Oliceridine is Associated with Reduced Risk of Vomiting and Need for Rescue Antiemetics Compared to Morphine: Exploratory Analysis from Two Phase 3 Randomized Placebo and Active Controlled Trials[END_REF]. At present, TRV130 is approved by the Food and Drug Administration (FDA) since august 2020 to treat moderate to severe acute pain. Finally, TRV130 turned to be a good analgesic medication but its advantages, especially regarding side effects, toward classical opioids such as morphine is still matter of debate.

Thanks to the progress in GPCR crystal structure, Manglik and co-workers used the MOR crystal structure to virtually dock millions of compounds and they discovered PZM21 [START_REF] Manglik | Structure-based discovery of opioid analgesics with reduced side effects[END_REF]. PZM21 was found to be biased towards G protein pathway and active in both spinal and supra-spinal analgesia in rodents [START_REF] Manglik | Structure-based discovery of opioid analgesics with reduced side effects[END_REF][START_REF] Kudla | Functional characterization of a novel opioid, PZM21, and its effects on the behavioural responses to morphine[END_REF].

Interestingly, this compound seemed to induce less respiratory depression and constipation than morphine [START_REF] Manglik | Structure-based discovery of opioid analgesics with reduced side effects[END_REF].

Finally, the latest notable MOR biased compounds are a series of a piperidine core structurecontaining molecules (SR series) generated from a synthesis-driven approach. Some of these compounds were found to be biased, at different levels, toward G protein. They induced analgesia (measured in the hot plate and tail flick tests), showed minor respiratory depression and are brain penetrant [START_REF] Schmid | Bias Factor and Therapeutic Window Correlate to Predict Safer Opioid Analgesics[END_REF]. Interestingly, SR-17018 was demonstrated to be more potent and efficacious than morphine or oxycodone in a chemotherapeutic-induced neuropathy pain model and displayed a weak tolerance after a repeated administration [START_REF] Pantouli | Comparison of morphine, oxycodone and the biased MOR agonist SR-17018 for tolerance and efficacy in mouse models of pain[END_REF]. Finally, it reversed morphine-induced withdrawal [START_REF] Grim | A G protein signaling-biased agonist at the μ-opioid receptor reverses morphine tolerance while preventing morphine withdrawal[END_REF].

But they remain MOR ligands (with their side effects)

Because they bind and activate MOR, all these biased ligands could still promote main effects attributed to MOR stimulation including reward system activation. Indeed, in rats, TRV130 displayed same effects as morphine in intracranial self-stimulation [START_REF] Altarifi | Effects of acute and repeated treatment with the biased mu opioid receptor agonist TRV130 (oliceridine) on measures of antinociception, gastrointestinal function, and abuse liability in rodents[END_REF]. In a drug-discrimination experiment conducted in rats, two biased ligands TRV130 and SR-14968 generalized to fentanyl [START_REF] Schwienteck | Effectiveness comparisons of G-protein biased and unbiased mu opioid receptor ligands in warm water tail-withdrawal and drug discrimination in male and female rats[END_REF], suggesting that these ligands may produce prototypic MOR agonist abuse-related effects. Using a model closer to the human behaviour, intravenous self-administration procedure in rats , Austin Zamarripa and colleagues demonstrated that TRV130 was equi-potent to oxycodone (Austin [START_REF] Zamarripa | The G-protein biased mu-opioid agonist, TRV130, produces reinforcing and antinociceptive effects that are comparable to oxycodone in rats[END_REF].

More recently, whereas PZM21 demonstrated a lack of effect in the conditioned place preference paradigm in two independent studies in mice [START_REF] Manglik | Structure-based discovery of opioid analgesics with reduced side effects[END_REF][START_REF] Kudla | Functional characterization of a novel opioid, PZM21, and its effects on the behavioural responses to morphine[END_REF], reinforcing effects, here again comparable to those induced by oxycodone, were demonstrated using a self-administration procedure in monkeys [START_REF] Ding | Antinociceptive, reinforcing, and pruritic effects of the G-protein signalling-biased mu opioid receptor agonist PZM21 in non-human primates[END_REF], but not in rats [START_REF] Kudla | Functional characterization of a novel opioid, PZM21, and its effects on the behavioural responses to morphine[END_REF]. All these data from animal studies could suggest that the MOR biased compounds could have reinforcing effects in humans. Indeed, in a clinical study assessing the analgesic effects of TRV130 in healthy volunteers, Soergel and collaborators evaluated the abuse-related subjective effects with a drug effects questionnaire. They found that the dose of 3 mg, equianalgesic to 10 mg morphine, induced the same subjective opioid effects such as "high" or "liking" [START_REF] Soergel | Biased agonism of the μ-opioid receptor by TRV130 increases analgesia and reduces on-target adverse effects versus morphine: A randomized, double-blind, placebo-controlled, crossover study in healthy volunteers[END_REF]. These data highlight the importance of the model at two level: the species and the experimental model used. Furthermore, in the same species, the choice of the strain of animals used can also have a very important impact on the results. Thus, the first results obtained by Bohn and co-workers on beta-arrestin 2 knockout animals were obtained on a mixed genetic background (beta-arrestin 2 knockout 129/SvJ backcrossed with C57BL6).

Previous studies reported that the 129/SvJ mice have a higher sensitivity to morphine-induced antinociception [START_REF] Mogil | Nociceptive and morphine antinociceptive sensitivity of 129 and C57BL/6 inbred mouse strains: Implications for transgenic knock-out studies[END_REF], did not develop tolerance to morphine-induced analgesia [START_REF] Crain | Enhanced analgesic potency and reduced tolerance of morphine in 129/SvEv mice: evidence for a deficiency in GM1 ganglioside-regulated excitatory opioid receptor functions[END_REF] and were less sensitive to respiratory depression [START_REF] Bubier | Genetic variation regulates opioid-induced respiratory depression in mice[END_REF]. When reassessing the pharmacological responses induced by opioids in betaarrestin 2 knockout using an homogenous C57BL6 background, Koblich and colleagues did not find any difference in tolerance to analgesic effects of few opioids including morphine [START_REF] Koblish | TRV0109101, a G Protein-Biased Agonist of the µ-Opioid Receptor, Does Not Promote Opioid-Induced Mechanical Allodynia following Chronic Administration[END_REF]. More recently, Kliewer and co-workers reported that the knockout of betaarrestin 2 in C57BL6 mice did not alter the effect of morphine on both respiratory depression and constipation compared to wild-type animals [START_REF] Kliewer | Morphineinduced respiratory depression is independent of β-arrestin2 signalling[END_REF]. Taken together, these data suggested that the use of 129/SvJ strain might have misled us on the role of beta-arrestin 2 in morphine effects and as a consequence it questions the G protein biased agonists as valuable strategy for safer analgesic medication. More importantly, an increase number of studies challenged the existence of the biased signalling itself (Gillis et al., 2020b). Evaluating the bias usually consist of measuring the ability of an agonist to stimulate two signalling pathways. At end, one will compare an amplified assay of G activation versus an unamplified assay of beta-arrestin recruitment. Amplified assay are relatively insensitive to agonist efficacy differences. Whereas an agonist with a low intrinsic efficacy will show a maximal effect in G protein assay, it will have a lower maximal effect in beta-arrestin assay. Using newly unamplified probes to measure signalling, several studies demonstrated that the newly biased ligands are in fact agonists with a low intrinsic efficacy relative to morphine (Gillis et al., 2020a;[START_REF] Stoeber | Agonist-selective recruitment of engineered protein probes and of GRK2 by opioid receptors in living cells[END_REF][START_REF] Vasudevan | Assessment of structure-activity relationships and biased agonism at the Mu opioid receptor of novel synthetic opioids using a novel, stable bio-assay platform[END_REF]. These recent data show that while biased ligands may have benefits over unbiased ligands, particularly in terms of the reduction of side effects, these may not be due to the biased nature of the compounds.

Concluding remarks

It took 30 years from the first observations of a biased signalling [START_REF] Jarpe | d-Arg1,d-Phe5,d-Trp7,9,Leu11]Substance P Acts as a Biased Agonist toward Neuropeptide and Chemokine Receptors[END_REF] and the release of Oliceridine, the first biased opioid receptor compound on the market. In the last few years, some studies have challenged this concept by suggesting that the observed bias is rather the expression of the agonist low intrinsic efficacies (Gillis et al., 2020a) or binding kinetics [START_REF] Herenbrink | The role of kinetic context in apparent biased agonism at GPCRs[END_REF]. Therefore, the promises of having opioid ligands with virtually a lack of side effects will still take time. Rather than a revolution, the biased compounds are a hopeful evolution in the pharmacological research of analgesics. The pharmacology of GPCRs is very complex, and the simplistic view of one receptor, one signalling pathway, one pharmacological response is a far too simplistic concept, which does not reflect the complexity of life sciences. So, should we abandon opiate ligands to treat pain? Certainly not, for now, as they remain the most effective antinociceptive agents. Nevertheless, alternative research strategies should be explored further. The development of non-selective opioid ligands (agonist or mixed agonist/antagonist) could be an interesting approach. For instance, cebranopadol, a mixed MOR/DOR/KOR/NOR agonist was found to be efficient in both acute and chronic pain with a delayed development of tolerance as compared to morphine [START_REF] Linz | Cebranopadol: a novel potent analgesic nociceptin/orphanin FQ peptide and opioid receptor agonist[END_REF]. This compound is actually in phase 3 trials. More recently, Ding and co-workers described the AT-121, a MOR/NOR mixed agonist with analgesic effects in non-human primates and a lack of the most frequent opioid-associated side effects such as physical dependence, abuse potential, respiratory depression, and opioid-induced hyperalgesia [START_REF] Ding | A bifunctional nociceptin and mu opioid receptor agonist is analgesic without opioid side effects in nonhuman primates[END_REF]. Allosteric ligands at the MOR could be also interesting as they could be used to enhance the effects of either the exogenous orthosteric agonists when administrated together, or the endogenous opioid peptides released in stress or pain condition [START_REF] Livingston | Allostery at opioid receptors: modulation with small molecule ligands[END_REF]. Another strategy could be the development of dualtarget ligands directed to both MOR and another receptor, including for instance dopamine D3 receptor [START_REF] Bonifazi | Novel Dual-Target μ-Opioid Receptor and Dopamine D3 Receptor Ligands as Potential Nonaddictive Pharmacotherapeutics for Pain Management[END_REF]. Whether biased ligands, ligands with low selectivity towards different opioid receptors, or bivalent ligands, the development of an analgesic

compound free of side effects and with low abuse potential seems to be a way still paved with some difficulties.

Author Contributions Statement

NM and FN wrote the manuscript. All authors take responsibility for final content. All authors read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Hind Sylvain Chabanet
email: sylvain.chabanet@univ-lorraine.fr

Bril El-Haouzi
email: hind.el-haouzi@univ-lorraine.fr

Philippe Thomas
email: philippe.thomas@univ-lorraine.fr

Toward a sawmill digital shadow based on coupled simulation and supervised learning models

Keywords: Digital shadow, Active Learning, Sawmill simulation

Digital Twins (DT) have been introduced as promising decision support tools in many different settings and to serve a variety of purposes. Many challenges are raised by their development, including an efficient usage of their computational resources to balance performance on precision, computational cost and speed. This study is, in particular, concerned with Digital Shadows (DS), a concept derived from DT, applied to sawmills sawing production units. A method to combine a computationally intensive sawmill simulation model with a machine learning model is proposed to predict set of lumbers sawed from logs. Numeric experiments are exposed, and the proposed method demonstrate improvements from 11% to 18% of the monitored couple regret from its baseline.

Introduction

The concept of digital shadows (DS) is derived from the concept of digital twin (DT) which is tracked back to NASA Apollo program, which built two identical vehicles, one sent to space missions and the other remaining on earth [START_REF] Boschert | Digital twin-the simulation aspect[END_REF]. The objective was to monitor the twin sent to space during mission and find solutions when faced with problems. Since then, the concepts and definitions proposed for DS and DT have evolved and it remains a field of active research. A broad generally accepted definition encompassing both DS and DT concepts was proposed as "a set of digital models of physical entities usable for optimization and decision-making purposes based on real-time information" [START_REF] Savolainen | Contrasting digital twin vision of manufacturing with the industrial reality[END_REF]. In particular, while both DS and DT are often associated with simulation, digital models based on many technologies have been propose to support their respective capabilities. For example, [START_REF] Gardner | Towards the development of an operational digital twin[END_REF] propose an hybrid DT based on physics-based simulation augmented by a machine learning (ML) model to predict acceleration responses of a structure subject to vibrations. Similarly, [START_REF] Magnanini | A model-based digital twin to support responsive manufacturing systems[END_REF] introduce a DT based on analytical models for the performance evaluation of a manufacturing system. While the difference between DS and DT is still subject to debates, an interesting discussion is provided by [START_REF] Sapel | Towards digital shadows for production planning and control in injection molding[END_REF]. DT should be able to automatically influence their physical counterpart. As such, they require an extremely precise representation and detailed models of its behavior. DS, on the other hand, are decision support tools designed to answer specific questions. It should be noticed, however, that this distinction is not always made in the scientific literature, and DS can hardly be studied independently from DT.

Due to the variety of technologies proposed to design them, DS might be better described based on their capabilities : to maintain a description of the physical counterpart (PC) based on real time or near real time data, to predict future states of the PC or compare various configurations, and to support decision making by returning instruction or optimized plans to the physical world. Such capabilities are naturally built on top of one another. Descriptive capabilities are required for predictive model to front-run production. Similarly, the output of predictive models can be fed as input to optimization models prescribing optimized plans (figure 1).

Fig. 1 The three main capabilities of a complete DS decision support tool Digital shadows and digital twins have been studied in very different contexts, from hospitals [START_REF] Karakra | Hospit'win: a predictive simulationbased digital twin for patients pathways in hospital[END_REF] to sawmills [START_REF] Chabanet | Toward digital twins for sawmill production planning and control: Benefits, opportunities, and challenges[END_REF]. Sawmills, in particular, might benefit greatly from DS for operational production monitoring and control. Sawing operations, the breaking of logs into lumber, is a process divergent and in co-production. From a single log, a sawmill will obtain simultaneously several pieces of lumber of varied dimensions (figure 2). In the remaining of this paper, the set of lumber obtained from a single log is named a basket of products (BoP). Many factors introduce uncertainty on the yield of the sawing process. These factors include, for example, the heterogeneity of the shapes of logs, whose general characteristics may change from one harvesting site to another, or the emergence of automated sawlines which adjust sawing operations on a log per log basis. All this makes it difficult to know in advance the mix of lumber that would be obtained from sawing a batch of logs. Solutions proposed in the literature to overcome this uncertainty include the use of specialized optimization methods, such as robust or stochastic optimization [START_REF] Zanjani | Sawmill production planning under uncertainty: Modelling and solution approaches[END_REF], or the usage of sawing simulators able to predict logs BoP based on information over their shape [START_REF] Maness | Multiple period combined optimization approach to forest production planning[END_REF] or even internal structure. Information over log shape may, in particular, take the form of 3D scans, for example obtained from laser scanners which may be placed at the entrance of sawmill log yards or sawing units. General simulators able to flexibly model many sawmill configurations and predict BoP of logs based on their scans have been developed by researchers and industrial. They include, for example, Sawsim1 or Optitek [START_REF] Goulet | Optitek: User's manual[END_REF]. Optitek, in particular, as been used in many academic studies [START_REF] Morin | Machine learning-based models of sawmills for better wood allocation planning[END_REF][START_REF] Wery | Simulation-optimisation based framework for sales and operations planning taking into account new products opportunities in a co-production context[END_REF]. Interestingly, both of these studies demonstrate the benefit of BoP predicted by either simulation or ML based metamodels of these simulators as input for planning decision support optimization models.

Such generic sawmill simulation tools may be great steps toward the development of flexible simulation-based sawmill DS. In particular, they would provide flexible yet precise models backing predictive capabilities. Several researchers, however, have noticed the important computational time associated with running these simulations for thousands of logs [START_REF] Morin | Machine learning-based models of sawmills for better wood allocation planning[END_REF][START_REF] Wery | Simulation-optimisation based framework for sales and operations planning taking into account new products opportunities in a co-production context[END_REF][START_REF] Morneau-Pereira | An optimization and simulation framework for integrated tactical planning of wood harvesting operations, wood allocation and lumber production[END_REF], which would impair their use for real or near real-time problems. The design of DS based on simulation and ML models used simultaneously, appears, therefore, as an interesting solution to this problem. Such DT and DS mixing data-driven methods and other simulation techniques are sometimes referred to as "hybrid" in the literature, for example by [START_REF] Asadi | Machine-learningenabled digital twin of welded structures for rapid weld sequence design[END_REF]. The "hybrid" qualifier refers only to the variety of models involved in the predictive capabilities of the twin. An interest of these hybrid DS is that, in particular, predictions performed by a ML model are, in general, very fast. Such models may require computationally intensive training phases, but these can be performed offline when necessary. [START_REF] Asadi | Machine-learningenabled digital twin of welded structures for rapid weld sequence design[END_REF], for example, develop such a model to combine a simulation model based on finite element analysis with a ML model to optimize welding operations.

Several factors complicate, indeed, the use of DS based only on data-driven, ML-based predictive models. ML models, in particular, require large datasets of labeled examples to be trained on. While in some cases these training examples may be fully obtained from the physical twin, this is complicated for sawmills by the difficulty of tracking logs and lumbers through several diverging transformation processes. In a more general setting, simulation remains necessary when historical data are not accessible from the start to train ML models, for example after changing machine settings. A second reason for favoring an hybrid DT is that ML models may, on some prediction problems and depending on the simulation technology, remain more prone to error and less explainable than these simulations.

This, however, raises the question of how to couple a simulation and ML models performing the same prediction task with their respective advantages and drawback. In this paper we propose a method to automatically decide, on a log per log basis, whether its BoP should be predicted by a ML model or a simulator.

The remaining of this study is structured as follows. Section 2 first presents in more details the hybridization problem considered. Section 3 details our proposed method, for which an experimental proof-of-concept is given in section 4. Section 5 concludes this paper and exposes future works.

Problem description

The general problem considered in this study assumes the existence of a stream of data items. For every individual data item X, some quantity y has to be predicted by a DS. In this study, a data item is a set of information over a specific log, and y is their BoP, but it could be generalized to any other type of stream. The existence of a DS containing at least two digital models able to make this prediction is similarly assumed. The first model is a simulation model, very precise but computationally intensive. The second is a ML model, fast but overall not as precise as a simulation. This ML model may be based on any supervised learning algorithm, for example neural network or random forest. As such, it is data-driven and requires labeled examples to be trained on. The quality of a model prediction ŷ is evaluated by a loss function l(ŷ, y), so that a prediction ŷ closest to the real value y will lead to a lower loss. No upper bound is assumed on the number of simulations that may be run in parallel. This may be an acceptable approximation if, for example, the simulation is run on a cloud service so that the amount of computational resources which can be punctually used is far greater than what is wished to be used on the long term to maintain overall computational costs reasonable. However, an asymptotic simulation budget b is targeted, i.e, a limit of the simulation budget as time goes to infinity. It is defined here as the limit of the total accumulated simulation time divided by the time length of the stream :

b = lim T →+∞ 1 T ∞ k=0 p k 1 {a k <T } . (1)
Here, p k is the time necessary to predict y for the k th data item entering the stream, a k its time of arrival and 1 {...} the indicator function.

The objective is to propose a dispatching function 1 T oSimul(X) taking value 1 if the prediction of the data item X has to be performed by the simulation model. Such a decision has to be taken immediately after arrival, without consideration for future data item. Additionally, the loss of a prediction from the ML model, l(ŷ M L , y) is, of course, unknown when the decision has to be taken and can't be used directly.

For the purpose of mathematical modeling, the following assumptions are also used. The prediction of the simulation model is supposed to be exact, i.e, l(ŷ Simul , y) = 0 if 1 T oSimul(X) = 1. The prediction time of the ML model is neglected, i.e, l(p k , y) = 0 if 1 T oSimul(X k ) = 0. The times of arrivals of data items are assumed to follow a Poisson process of rate µ. Simulation times are sampled from an Exponential law of mean λ. These choices make the process recording the number of simulations being run simultaneously a continuous-time markov chain process, ensure the existence of a stationary limiting distribution and allows the usage of the ergodic theorem. However, the results described in the following of this section hold in a more general case, as long as stationarity and ergodicity can be assumed.

Following this modeling, consider α = P(1 T oSimul(X) = 1) the fraction of the stream whose prediction is obtained from the simulation model. If the decision 1 T oSimul(X k ) is taken independently from past and future arrival, the process formed by the data items for which a prediction is requested from the simulation model forms a Poisson process of rate µα. Noting N (t) the number of simulations running in parallel at time t, the budget b can be rewritten as:

b = lim T →+∞ 1 T T 0 N (t) dt , (2)
i.e., the target budget b is also the average number of simulation running in parallel at any time. Additionally, according to the ergodic theorem, it is equal to E(N ) under the process stationary distribution. Little's law [START_REF] Little | Little's law[END_REF] yield, in turn:

E(N ) = E(p)µα = λµα . (3)
Hence, to target the budget b, the simulation rate α must be α = b λµ .

Proposal

The problem considered in this study is similar to problems detailed in the active learning (AL) literature. AL is a field from the ML literature concerned with reducing the cost associated with labeling data to train an ML model. While unlabeled collection of data items X are, in many settings, easy to gather in huge quantity, their labels y may be costly or difficult to obtain. This is often the case when the Oracle providing the labels is, for example, a team of human experts. More precisely, AL is concerned with guiding the selection of a subset of unlabeled items that are considered the more useful to improve the performances of the ML model. It is often an iterative process where a ML model has already be initially trained on a small training set and is used to define a measure m of the utility of labeling a point. Data items maximizing this measure, or having it above some threshold are then selected from an unlabeled dataset or data stream and labeled. The ML model may be periodically retrained when a sufficient amount of new labeled data have been gathered, or when another condition is trigered such as the detection of changes in the data stream. For example, [START_REF] Weigl | On improving performance of surface inspection systems by online active learning and flexible classifier updates[END_REF] propose an AL method based on a Random Forest (RF) model to detect surface defects in a stream of products.

Three main scenarios exist in the AL literature, depending on the way the active learner gains access to data :

• Pool-based active learning is concerned with cases where an unlabeled dataset may be accessed at once by the learner. • Stream-based active learning considers cases where data items are not accessible at once but one by one from a data stream. The decision to label every item has then to be taken immediately for every data item, or for small batches of item, without considering the future of the stream. • Membership Query Synthesis scenarios do not consider the existence of real unlabeled data. Instead, the learner may create synthetic data and request labels for these. This scenario is rarely investigated as synthetic data may be difficult to label, especially for human oracles.

All three scenarios may be of interest in the context of DS. In particular, simulation-based DS may be of great interest to label synthetic data and stream-based active learning may be used to manage real time data streams continuously collected by the twin. This study focuses on this second stream based case.

Many measures m of labeling utility have been introduced for both regression and classification tasks [START_REF] Kumar | Active learning query strategies for classification, regression, and clustering: a survey[END_REF]. In particular several are based on heuristic or statistical measures of the uncertainty or confidence associated with a prediction. The method proposed in this study is to use such measures of uncertainty to design the dispatching function 1 T oSimul(X) , under the reasoning that such a measure may bring valuable insight over the unobservable loss l(ŷ M L , y), and that ML prediction whose confidence is considered too low should preferably be obtained from the simulation model.

In particular, consider q the quantile at level 1 -α of m(X). In general, this quantile may be estimated from observed past values from the stream. If the dispatching function is taken as 1 T oSimul(X) = 1 m(X)>q , then a simulation result is requested for the portion α of the stream for which the ML model is considered the least confident according to m.

The ML model used in this study is a RF model [START_REF] Breiman | Random forests[END_REF]. RF are machine learning models popular for their ease of use, overall good results and robustness, and have been used with success in many studies in industrial context [START_REF] Ruiz | Machine learning algorithms for the prediction of the strength of steel rods: an example of datadriven manufacturing in steelmaking[END_REF][START_REF] Soto | An online machine learning framework for early detection of product failures in an industry 4.0 context[END_REF]. They also gave good performances in past studies comparing ML models to predict BoP of logs [START_REF] Morin | Machine learning-based models of sawmills for better wood allocation planning[END_REF][START_REF] Morin | Machine learning-based metamodels for sawing simulation[END_REF]. RF are ensembles of decision tree models, where the prediction is given by ŷ = t(X) = 1 B B i=0 t i (X). Here, t i (X) are predictions of individual trees and B is the number of trees in the forest. Individual trees are trained on bootstrap samples selected from the whole training data set.

The measure m considered here is an estimator of the sampling variance of a random forest introduced by [START_REF] Wager | Confidence intervals for random forests: The jackknife and the infinitesimal jackknife[END_REF]. It evaluates how much the prediction may change if a similar model had been trained on a different dataset. This measure was used for AL for example in [START_REF] Farache | Active learning and molecular dynamics simulations to find high melting temperature alloys[END_REF] which propose a method to predict alloys melting points. It is expressed as :

k i=1 Cov[J i , t . (x)] 2 + k i=1 ( t-i (x) -t(x)) 2 - ek B 2 B i=1 (t i (x) -t(x)) 2 . ( 4
)
Here, k is the member of data items in the training set, e the Euler number, t-i (x) the average of the prediction of trees which didn't use the i th data item for training, and J i the number of time the i th item was present in the bootstrap sample used to train a tree t . . Cov is the usual covariance estimate. Considering that [START_REF] Wager | Confidence intervals for random forests: The jackknife and the infinitesimal jackknife[END_REF] recommend to use forests with high number a trees, comparable with the number of data items used to train it, B is fixed to 500 in this paper experiments.

If the simulation model is precise enough, its predictions, collected along the stream may, additionally, be added to a training dataset used to periodically retrain and improve the performances of the ML model with more data.

A pseudo-code summary of the proposed solution is exposed algorithm 1.

Algorithm 1 General flow of the proposed strategy

Experiments and results

This section first presents the dataset used during the experiments proposed in this study. Additional details over these experiments are also given, and the results are exposed.

Materials and method

The dataset used in the experiments described in this section originates from the Canadian forest-product industry. It contains information about 2219 wood log. In particular, both the 3D scans of every logs given, as well as a collection of six describing features. An example of such a 3D scan is presented figure 3. A scan point cloud, whose points are organized into rough ellipsoids spanning the log surface. The six know-how features summarizing information over the shape of logs are the length, volume, diameter at both extremities, curvature and taper (a measure of the shrinking of a log from one end to another). While the scans were used to simulate the sawing of logs and obtain their BoP, these six knowhow are required by the ML model to make their own predictions. 3D scans are, indeed, unordered collection of points and the number of points varies from scan to scan. They are, therefore, unstructured data and can't be used directly by commons ML models. To summarize scans by a few descriptive features implies a loss of information that cannot be used by the ML model to make predictions. Therefore, such a ML model cannot, in particular, be expected to perform as well as a simulation model.

The BoP of every log, obtained by simulating its sawing with the simulator Optitek were also obtained. Simulation times, however, were not part of this dataset. They have been, therefore, generated randomly when required.

The sawmill simulated was able to produce up to 47 types of lumbers, characterized by their length, width and thickness. For this reason, BoP are represented as vectors of dimension 47, where the i th element encodes how many pieces of lumber of type i are present in the BoP. The machine learning prediction task is, therefore, modeled as a multi-output regression problem. Considering that the measure described equation 4 is initially designed for single output regression problems, in practice a collection of 47 measures m i where used, one for every type of lumber, and their sum m(X) = 47 i=1 m i (X) was used.

To simulate a stream for experimental purpose, the dataset is ordered at random. The first 500 elements are used to initialize the RF model. Interarrival time in the stream were sampled from an exponential distribution with average λ = 1 and simulation time were sampled from an exponential distribution with rate µ = 1 10 so that the average simulation time is in average 10 time greater than the average inter-arrival time. Several values of the budget were tested.

At the end of the stream, the total stream time and accumulated simulation time are collected to estimate the effective simulation budget. The fraction of the stream sent to the simulation model is also collected. Two evaluation scores are additionally collected. The first one is the ML regret:

1 n n l=1 l(ŷ M L,l , y l ) (5)
with n the number of data items in the stream, y l the label of the l th arriving item, and ŷML,l its prediction performed by the ML model whether it is sent to the simulator or not. This score evaluates the performance of the ML model if it were used for every prediction.

Similarly, the coupling regret is defined as:

1 n n l=1 l(ŷ M L,l , y l )(1 -1 {T oSimul(X l )} ) (6)
and measures the performance of the hybridization. The performances of the proposed dispatching function will, additionally, be compared with the performances of a random dispatching requesting simulation depending on the result of a collection of i.i.d Bernoulli variables of parameter α. The ML model was retrained every time 100 new labeled items, or, in average, every 400 new arrivals.

To average out the impact of the exact ordering of a stream, and time sampling, this will be repeated 50 times with different random seeds.

Results

Results for three different values of the budget 1, 2.5 and 4 are presented in table 1. In particular, this means that, according to section 3, the use of the simulation model should be requested for approximately 10%, 25% and 40% of the logs from the stream to predict their BoP. This is, indeed what is observed in practice and is recorded in the column "simulated ratio". Additionally, the respective targeted asymptotic budgets are respected, as can be seen from the observed accumulated simulation times. It is interesting to notice that the ML regret, which evaluates the performance of the ML model if it were used for all predictions, is always slightly lower for the AL-based method than for the random method. It is, indeed, the initial objective of AL methods to improve the performances of ML models by guiding the selection of data items that should be labeled rather than sampling at random. Performing Welsh t-test to compare the averaged ML regrets over the various repetition of the experiment yields p-values of 0.001 for b = 1, 1.4 × 10 -7 for b = 2.5 and 1.8 × 10 -7 for b = 4. In every case, the null hypothesis of equality of the samples means can, therefore, be rejected. Despite being small, this difference is, therefore, considered statistically significant, w.r.t the repetition of the experiment. This includes randomness of the stream order, simulation time and inter-arrival time. It should be notice that the randomness from the original data collection process itself is not considered by the test. These ML regrets remain, however, high with respect to the average number of products in a basket, at 4.6, and the quadratic error obtained by a predictor predicting always the average BoP over the dataset, at 8.3. This might be partly caused by the small size of the datasets used to train these predictors.

Coupling regrets are, of course, always lower than ML regret, due to a fraction of the prediction being requested from the simulator and, therefore, considered to be exact. Once again, however, the coupling regret of the ALbased method is lower than the one of the random method, to a greater extent than what was already observed in the case of the ML regret. In particular, the regret is reduced by 11% for b = 1, 15% for b = 2.5 and 18% b = 4. Therefore, the sampling variance measure used in the AL-base method appears efficient at anticipating predictions from the ML model with higher loss, and requesting predictions from the simulation model instead.

Conclusion

Many challenges remain in the development of DS in general, and sawmill DS in particular. The development of methods to allow an efficient usage of computational resources balancing computational cost with analytic performances appears, in particular, important for both economic and sustainability reason. To this end, the development of DS based on hybrid simulation-ML models able to benefit from the advantages of both technologies while compensating their drawback.

This study introduces a method, based on AL, to couple a computationally intensive sawmill simulator with a less precise but fast ML model to predict BoP of wood logs. The main objective is to make the best use of the computational resources allocated to the simulation model by requested simulations only for these logs for which the ML model predictions appear the least certain according to some statistical measure. During experiments, this proposed method demonstrated reductions of the simulation ML couple regret ranging from 11% to 18% depending on the simulation budget. Additionally, to use the simulation results to improve the performances of the ML model itself demonstrate a slightly better effect with the proposed AL-based method than with a random dispatching of the prediction toward the simulator.

However, several limits exist to the results presented in this paper. Firstly, more experiments need to be run to assess the impact of various parameters, in particular the average simulation time and stream arrival rate. Secondly, the respect of the simulation budget is asymptotic and based on statistical assumptions, including stationary of the data stream. This mean, in particular, that this budget may be temporarily exceeded, and increase significantly if drifts were to happen in the stream. Therefore modifications of the proposed strategy enforcing the budget to always remain between acceptable bounds appear particularly important. Additionally, many more utility measures have been proposed in the AL literature. Further works comparing different such methods will therefore be necessary.

Fig. 2

 2 Fig. 2 Sawing a log produce several pieces of lumbers with various dimensions

Fig. 3

 3 Fig. 3 3D scan of a log

Table 1

 1 Average and standard error over the 50 experiments of the simulation time over stream time, ratio of the stream sent to the simulator, and regrets

	budget Method	simulation time ratio	simulated ratio	ML regret	coupling regret
	1	Random method	1.01 (0.10) 0.10 (0.01) 3.25 (0.13) 2.92 (0.12)
		AL-based method	1.12 (0.07) 0.11 (0.004) 3.17 (0.10) 2.59 (0.08)
	2.5	Random method	2.49 (0.17) 0.25 (0.01) 3.16 (0.09) 2.38 (0.09)
		AL-based method	2.46 (0.15) 0.25 (0.004) 3.07 (0.07) 2.03 (0.06)
	4	Random method	4.00 (0.21) 0.40 (0.01) 3.09 (0.08) 1.85 (0.07)
		AL-based method	3.96 (0.17) 0.39 (0.01) 3.01 (0.06) 1.52 (0.05)

https://www.halcosoftware.com/software-1-sawsim, Last accessed on June,

Acknowledgment

The authors gratefully acknowledge the financial support of the ANR-20-THIA-0010-01 Projet LOR-AI (lorraine intellgence artificielle) and région Grand EST. We are also extremely grateful to FPInnovation who gathered and processed the dataset we are working with.