Robin Tournemenne
email: robin.tournemenne@inria.fr

Jean-Franc ¸ois Petiot
email: jean-francois.petiot@ls2n.fr

Bastien Talgorn

Jo Ël Gilbert
email: joel.gilbert@univ-lemans.fr

Michael Kokkolaras
email: michael.kokkolaras@mcgill.ca

Sound simulation-based design optimization of brass wind instruments

Brass Instrument optimization

We present a method for optimizing the inner shape of brass instruments using sound simulations. This study considers different objective functions and constraints (representative of both the intonation and the spectrum of the instrument) for a relatively large number of design variables. A complete physics-based model, taking into account the instrument and the musician embouchure, is used to simulate permanent regimes of sounds by means of the harmonic balance technique, the instrument being represented by its input impedance. The design optimization variables are related to the geometrical dimensions of the resonator. The embouchure's parameters are varied during the optimization procedure to obtain an average behavior of the instrument. The objective and constraint functions of the optimization problem are evaluated using the physics-based simulation model, which is computationally expensive. Moreover, the gradients of the objective and constraint functions can be discontinuous, unavailable, or hard to approximate reliably. Therefore, we employ a surrogate-assisted derivative-free optimization strategy using the mesh adaptive direct search algorithm (MADS). One example of a B trumpet's bore is used to demonstrate the effectiveness of the design optimization approach: the obtained results improve previously reported objective function values significantly.

I. INTRODUCTION

The development of innovative and higher-quality designs is crucial to the viability of musical instrument manufacturers. Many prototypes are required to ensure that quality attributes such as intonation, ease of emission, timbre, projection, etc. are adequate at each stage of the development process. Numerical acoustics models may be helpful in shortening development cycles by minimizing the manufacturing of costly prototypes. For example, in the case of brass instruments, several studies propose the use of numerical modeling to predict quality and guide the design process [START_REF] Campbell | Brass instruments as we know them today[END_REF][START_REF] Macaluso | Trumpet with near-perfect harmonicity: design and acoustic results[END_REF].

The dominant physical quantity impacting the sound quality of a brass instrument is its input impedance (cf. Figure 1). Input impedance is the frequency-dependent quotient of pressure and of the instrument volume flow at the instrument entry plane, and, at first approximation, is the result of the in-strument's interior shape (the bore). Many works have focused on modeling input impedance, e.g., [START_REF] Caussé | Input impedance of brass musical instruments-comparison between experiment and numerical models[END_REF].

The impact of input impedance on sound quality is acknowledged in [START_REF] Campbell | Brass instruments as we know them today[END_REF], where, in a first approximation attempt, playing frequencies are governed mainly by corresponding impedance peaks [START_REF] Eveno | The relationship between bore resonance frequencies and playing frequencies in trumpets[END_REF]. Beyond influencing the instrument intonation, impedance also determines an instrument's timbre and playability, notably due to the height and bandwidth of peaks.

With this in view, several researchers have used input impedance to design an instrument's inner shape using an optimization approach. Following Kausel's successful reconstruction of a trumpet bore using the Rosenbrock algorithm [START_REF] Kausel | Optimization of brasswind instruments and its application in bore reconstruction[END_REF], Noreland optimized the instrument's intonation with a hybrid scheme for the input impedance model and shape constraints [START_REF] Noreland | A hybrid scheme for bore design optimization of a brass instrument[END_REF]. Braden also optimized the intonation and the input impedance peak heights of a trombone using a multi-modal input impedance model [START_REF] Braden | Trombone bore optimization based on input impedance targets[END_REF], while Macaluso optimized and built a near-perfect harmonic trumpet [START_REF] Macaluso | Trumpet with near-perfect harmonicity: design and acoustic results[END_REF]. Some studies investigated the relationship between input impedance features and psycho-acoustic criteria: Poirson optimized the trumpet using objective functions based on the input impedance and targets defined by trumpet players preferences [START_REF] Poirson | Integration of user perceptions in the design process: Application to musical instrument optimization[END_REF]. Guillauteau looked for empirical relations between playing frequencies and resonance frequencies to optimize clarinets [START_REF] Guilloteau | Conception d'une clarinette logique[END_REF].

However valuable, these works focused exclusively on instrument performance, neglecting a crucial element in sound production: the musician embouchure. In particular, the studies of Eveno et al. showed that the relationship between the resonance frequencies of the impedance and the actual frequencies of the sounds played by musicians can vary significantly [START_REF] Eveno | The relationship between bore resonance frequencies and playing frequencies in trumpets[END_REF].

Although the impedance of an instrument provides interesting information about sound quality, prediction of "playability" and sound qualities of brasses based solely on impedance remains difficult.

A second approach is based on a holistic model of the physical phenomenon, coupling the instrument and the musician embouchure, to produce sound simulations representative of the instrument quality. Using this approach, the authors integrate sound simulations in the optimization process. These simulations are obtained from a physics-based model to account for the interaction of the instrument with a virtual musician embouchure [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF]. In a previous study, the instrument's intonation was optimized based on simulated playing frequencies [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF]. Two examples were considered, optimizing 2 and 5 of the bore's geometrical parameters, respectively; results were quite encouraging.

The main objective of the present paper is to extend this new optimization paradigm in order to assess both its potential and limitation. The two main novelties of this paper are the optimization of criteria based on the instrument sounds spectra, and the inclusion of constraints in the problem formulation. Three other contributions are noticeable: i) an improved version of the optimization method has been implemented, ii) a new solver is introduced, and iii) the performance limits of the optimization method are tested by considering 10 geometrical variables of the bore. A trumpet is used as a representative brass instrument to demonstrate the proposed design optimization method.

The paper is organized as follows. We first present extensive details on the physics-based model and the simulation technique. We then formulate the optimization problems and describe the principles of the MADS algorithm and the framework for surrogate-assisted optimization.

Finally, we conduct a case study concerning the shape optimization of a trumpet with ten design variables and draw conclusions.

II. TRUMPET MODELING

In this study, we utilize an elementary model of a brass instrument under playing conditions:

the vibrating lips are modeled as a one-degree-of-freedom (1-DOF) outward-striking valve, nonlinearly coupled to the air column of the brass instrument. This elementary model is a good compromise between simplicity and efficiency. While the 1-DOF model cannot model real musician lips exactly, (Yoshikawa, 1995) it is able to mimic a large range of playing phenomena (see for example the pioneering work of Elliot and Bowsher (1982), [START_REF] Elliott | Regeneration in brass wind instruments[END_REF] or the more recent works of [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF] [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF] and [START_REF] Velut | How well can linear stability analysis predict the behaviour of an Outward-Striking valve brass instrument model?[END_REF] [START_REF] Velut | How well can linear stability analysis predict the behaviour of an Outward-Striking valve brass instrument model?[END_REF]). Similar to [START_REF] Chen | Nature of the lip reed[END_REF], [START_REF] Chen | Nature of the lip reed[END_REF] we argue that while the lips may not be entirely modeled by a 1-DOF model, most characteristic behaviors of brasses can be reproduced by an outward striking reed.

Our physics-based model of the trumpet is based on Equations (1), (2), and (3), which all depend on three periodic variables: the opening height h(t) of the two lips, the volume flow u(t) of the air jet through the lip channel and the pressure p(t) in the mouthpiece (cf. Figure 2). of the opening height of the lips) and Q (the quality factor of the resonance of the lips). Input impedance Z is computed using the transfer matrix method considering plane wave propagation, visco-thermal losses [START_REF] Caussé | Input impedance of brass musical instruments-comparison between experiment and numerical models[END_REF], and a radiation function under an infinite plane baffle hypothesis.

p( jω) = Z( jω)û( jω) (1)
d 2 h(t) dt 2 + 2π f Q dh(t) dt + (2π f ) 2 (h(t) -h 0 ) = P m -p(t) µ (2) u(t) = bh + (t) sign(P m -p(t)) 2|P m -p(t)| ρ (
It is important to assess the validity of this 1-DOF model by indicating which behaviors of brasses can be reproduced adequately with simulations and which ones cannot. Previous results using this elementary model [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF] showed that the trumpet sounds simulated are dissimilar to the real trumpet sounds played by a musician. In particular, sound spectra in permanent regime are very different. This can be explained by inherent limitations of the model. A first limitation concerns the linear approximation of sound propagation in the brass resonator defined by its acoustic impedance (Equation ( 1)). When the instrument is played loudly and brassy, this approximation is no longer valid, and the nonlinear propagation needs to be taken into account. [START_REF] Myers | Effects of nonlinear sound propagation on the characteristic timbres of brass instruments[END_REF] A second limitation concerns the lip model, which simplifies the complicated real-lips motion [START_REF] Bromage | Open areas of vibrating lips in trombone playing[END_REF][START_REF] Martin | Directivity and the acoustic spectra of brass wind instruments[END_REF]Yoshikawa, 1995) with a 1-DOF (Equation ( 2)). Two-degrees-of-freedom (2-DOF) models have been considered for time domain simulation [START_REF] Adachi | Trumpet sound simulation using a two-dimensional lip vibration model[END_REF][START_REF] Boutin | Relationships between pressure, flow, lip motion, and upstream and downstream impedances for the trombone[END_REF]. Furthermore, measured mechanical responses of artificial [START_REF] Cullen | Brass instruments: Linear stability analysis and experiments with an artificial mouth[END_REF] and real brass-players lips [START_REF] Newton | Mechanical response measurements of real and artificial brass players lips[END_REF] revealed that a pair of mechanical resonances requires a 2-DOF model in order to be consistent with near threshold oscillations [START_REF] Cullen | Brass instruments: Linear stability analysis and experiments with an artificial mouth[END_REF]. However, although additional terms can theoretically be added to the 1-DOF model, the difficulty lies in selecting realistic values for the additional parameters [START_REF] Velut | How well can linear stability analysis predict the behaviour of an Outward-Striking valve brass instrument model?[END_REF].

A third limitation relates to the assumption that the volume flow u (Equation (3)), which controls the valve effect, is proportional to the opening height h between the lips. Experimental data reported in [START_REF] Bromage | Open areas of vibrating lips in trombone playing[END_REF] for a large set of playing frequencies and sound levels showed that the relationship may be exponential instead of linear, and dependent on the pitch and dynamic level of the note played.

Nevertheless, even if nonrealistic for the spectrum of trumpet sounds, previous studies confirm that this elementary model behaves in agreement with the main physical principles that govern the playing of brasses [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF][START_REF] Poirson | Study of the brightness of trumpet tones[END_REF]. In particular, results presented in [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF] showed that the elementary model is able to produce differences between instruments according to playing frequency, spectral centroid, and evolution of the spectral centroid with the playing dynamics that are, on average, in agreement with the differences noticed when a real trumpet is playing. This results justify the use of this elementary model in an optimization process for objective functions based on intonation, spectral centroid, or the evolution of the spectral centroid.

Numerical solutions of this system of equations are obtained using the harmonic balance technique to simulate the sound created by a given trumpet (defined by its input impedance Z) for a given "virtual musician embouchure" (defined by its control parameters). The harmonic balance technique considers the sound's permanent regime (steady state); since the latter is periodic, the truncated pressure is given by

p(t) = A 0 + N n=1 A n e i2πnFt + A * n e -i2πnFt . (4)
The unknowns, i.e., the amplitudes of the harmonics A n and the playing frequency F, are determined using Newton's method [START_REF] Gilbert | Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique[END_REF].

To perform a sound simulation, it is necessary to define relevant values (i.e., values that lead to a convergence towards a steady-state sound for a given note) for the parameters of the musician embouchure. For a given note, experience shows that countless embouchures may lead to a steadystate note. The choice of the range of the parameters is based both on numerical tests of the simulations and on measurements of real trumpet players. The three variables P m , µ , and f are considered as control parameters of the simulations, and constitute the virtual embouchure. The pressure P m in the mouth influences mainly the dynamics of a simulated sound and ranges from 1 to 12 kPa [START_REF] Fletcher | Blowing pressure, power, and spectrum in trumpet playing[END_REF]). In the following numerical experiments, we partitioned this range into three parts running from 1 to 5 kPa for what we call piano (p) dynamics, 5 to 9 kPa for mezzoforte (m f ) dynamics and 9 to 12 kPa for fortissimo ( f f ) dynamics. The frequency of the lips f enables the selection of the played regime (note): the higher the value of f , the higher the simulated regime. Exploration tests led to a range for f that spans from 130 Hz to 480 Hz to simulate the 2 nd , 3 rd , 4 th and 5 th regime of the B trumpet with no valve pressed, the regimes considered in this study. These regimes correspond to the musical notes B 3, F4, B 4, D5-concert-pitch. Finally, in order to produce many different sounds for every regime, we add variability to the embouchure making µ a control parameter of the simulations, ranging from 1 to 6 kg/m 2 [START_REF] Cullen | Brass instruments: Linear stability analysis and experiments with an artificial mouth[END_REF]. In our study, the values of b, Q , and h 0 are the same for every simulation [START_REF] Cullen | Brass instruments: Linear stability analysis and experiments with an artificial mouth[END_REF].

The values of the control parameters considered in this study are summarized in Table I.

Given that above 3000Hz the impedance magnitude is flat (see Figure 1), it is not relevant to consider many harmonics for the sound simulation. The highest studied note being D5 ( 587Hz), we chose to simulate our permanent regime with only N = 6 harmonics. In conclusion, for a given trumpet (characterized by its input impedance Z) and for a virtual musician embouchure (characterized by the parameters P m , µ , f , b, h 0 , and Q ), the simulation may generate one note (if the system converges), corresponding to one of the regimes 2, 3, 4 and 5 of the trumpet. Each note is characterized by its playing frequency F and the complex amplitudes of its 6 first harmonics.

It is important to mention that the computed sound p(t) corresponds to the sound in the mouthpiece. According to Benade [START_REF] Benade | Relation of AirColumn resonances to sound spectra produced by wind instruments[END_REF], a relevant spectrum transformation function could be defined to compute the sound outside the instrument. The difficulty lies in the definition of the radiated pressure, relevant from a perceptual point of view. For the optimization considered in this work, the well-defined pressure in the mouthpiece is deemed sufficient. It is also important to mention that the convergence of the simulation toward auto-oscillations is not ensured for a given shape of the resonator and embouchure. The search of convenient embouchures is a complex task, described in the following section.

III. OPTIMIZATION PROBLEM FORMULATION

The design optimization problem of an instrument can be formulated as the search for the optimal geometry minimizing an objective function:

min x∈Ω J(x), (5)
where J : R n → R is the objective function, and the vector x ∈ R n includes the design optimization variables. The design space Ω is a subset of R n delimited by box (bound) constraints. The design optimization variables are the geometric parameters that define the inner shape of the bore. To facilitate the input impedance calculations, the bore is approximated by a series of conical and cylindrical waveguide segments. Consequently, x is a vector of geometric quantities such as the lengths and radii of cylinders or cones. The design space Ω may be modified to obtain viable trumpet shapes.

Two classes of objective functions are available considering our physics-based model and the harmonic balance technique: descriptors based on playing frequencies and descriptors based on sounds spectra. Many quantities based on frequency and sound spectrum can be found in the literature; however, there is no consensus in the community regarding their influence on the instrument's musical quality. These disagreements notwithstanding, this paper considers three different objective functions based on intonation and spectral centroid given their recognized impact on the instrument quality [START_REF] Deutsch | The Psychology of Music[END_REF]. (see Section IV for details). The new design candidate x is evaluated in the same manner, and the process is iterated until the optimization algorithm termination criterion is satisfied.

Objective function based on intonation

For each note i, the average (across different embouchures) playing frequency F i (x, ϕ) at mezzoforte dynamics is computed. The intonation of the note is assessed by the deviation of the actual playing frequency, as simulated using the physics-based model, from the expected playing frequency. This objective function relies on a reference note from which the ideal and actual musical distance is computed. We chose for reference the 4 th regime of the trumpet with no valve pressed (B 4, concert pitch), given that it is the usual tuning note of the instrument. The equal-tempered scale is used to define the ideal distance between the studied note i and the reference note (Tournemenne et al., 2017).

For every note i, we compute the equal-tempered deviation (ETD) between the average frequency of the i th note F(x, ϕ i ) and the reference frequency F(x, ϕ r ) as

ETD i (x, ϕ) = α r→i -1200 log 2       F(x, ϕ i ) F(x, ϕ r )       , (7)
where α r→i is the ideal difference between the reference note r and the targeted note i given by the equal-tempered scale (-500 cents for example between B 4 and F4). The objective function J 1 (x, ϕ) for the whole instrument is the average of the absolute deviation across the (P-1) notes (notice that the deviation between the reference note and the 4 th note is always equal to zero, given that it is the tuning note)

J 1 (x, ϕ) = 1 P -1 i∈Notes ETD i (x, ϕ) . ( 8
)
2. Objective function based on the spectral centroid

The average spectral centroid SC(x, ϕ) at mezzoforte dynamics for every note is computed according to the 6 harmonic amplitudes A n of each simulated sound

SC(x, ϕ) = 1 card(ϕ) E v ∈ϕ SC(x, E v (P m , µ l , f l )) , (9)
where E v is an embouchure of the virtual musician with

SC(x, E v ) = 6 n=1 |n A n (x, E v )| 6 n=1 |A n (x, E v )| . (10)
Consequently, SC(x, E v ) spans from 1 to 6, representing the normalized spectral centroid. In this work, we decided to look for the instruments having the highest spectral centroid. These instruments would generally be considered as bright by musicians [START_REF] Poirson | Study of the brightness of trumpet tones[END_REF]. This is a somewhat arbitrary choice and other relevant descriptor/target may be found, although such consideration is out of the scope of the paper.

Consequently, the objective function J 2 (x, ϕ) is

J 2 (x, ϕ) = SC(x, ϕ). ( 11
)

Objective function based on the spectral centroid dynamics

This descriptor represents the ability of the instrument to maximise the spectral centroid difference between a piano (p) and a fortissimo ( f f ) dynamics. The idea is to find the instrument producing bright notes for high dynamics (high SC) while keeping a dark sound (low SC) for low dynamics. The average spectral centroid for each note and each piano (SC p (x, ϕ i )) and fortissimo (SC f f (x, ϕ i )) dynamics are computed. The objective function is then formulated as

J 3 (x, ϕ) = i∈Notes SC f f (x, ϕ i ) -SC p (x, ϕ i ) P = i∈Notes ∆SC i P . ( 12
)
For this descriptor, the note D5 concert-pitch has been discarded because of the difficulty in simulating it for low dynamics. A more application-oriented study considering the entire instrument's tessitura should account for these kind of difficulties.

B. Optimization problems

In summary, we consider three different optimization problems labelled Int, SC Int, and SC Dyn. They represent instrument design problems that are realistic from a musician point of view:

• Int: intonation improvement:

min x∈Ω J 1 (x) (13)
• SC Int: spectral centroid improvement under intonation constraint:

max x∈Ω J 2 (x) subject to J 1 (x) ≤ J 1 max (14)
• SC Dyn: spectral centroid dynamics improvement:

max x∈Ω J 3 (x) (15)
The constraint function of the SC Int problem is the objective function of the Int problem. It represents the scenario where musicians are willing to trade some of their instrument's intonation quality for a more brighter timbre.

The main challenge lies in the simulation of many different sounds (represented by their permanent regimes) during the numerical evaluation of the objective function. Practically, for each geometry x, we need to find suitable virtual embouchures leading to convergence toward permanent regime. Furthermore, similar to an inexperienced player that would blow even the most well-designed trumpet with a terrible sound, the virtual embouchure must be carefully selected in order to produce realistic sounds. No analytical approach exists to deal with this challenge and simple solutions always fall short. For example, it is not possible to define, a priori, a fixed list of virtual embouchures that will be used for every geometry x, because experience shows that the intersection of the sets of virtual embouchures leading to convergence toward a permanent regime for each geometry may be empty. It is far too expensive to process a complete fine grid of the 3 virtual embouchure parameters for every geometry x. Consequently, a rigorous preprocessing of the simulations is undertaken to help the simulations obtain a set of appropriate embouchures that converge toward adequate sounds for every geometry x.

This preprocessing is based on an exploration of the area of the design space Ω augmented by the 3 embouchure parameters leading to convergence of the sound simulation. If x is in R 2 , the space to explore has 5 dimensions: 2 geometric variables and 3 embouchure variables (P m , µ l , f l ). To explore this space, a five-dimensional Latin hypercube is built and the harmonic balance technique tries to simulate every sample.

In order to discard the simulations of unrealistic sounds mentioned above, we use a criterion representing the amplitude of the simulated sound relatively to the pressure in the mouth. If the amplitudes of the harmonics are large enough relative to the mouth pressure P m produced by the virtual musician, the sound is considered appropriate

6 n=1 A 2 n P m ≥ α. ( 16
)
Given the exploration and the criterion, an empirical method finds adequate virtual embouchures for any geometry x of the design space. There are two key differences between the preprocessing presented here and that reported in [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF]: i) the definition of the criterion threshold α based on live recordings and ii) a more robust technique used to define the set of virtual embouchures ϕ for any bore, both summed up in the following paragraph.

The live recordings of 3 helped trumpeters playing several times the 4 notes allowed us to estimate the standard deviation of playing frequency for each note . We then defined α in order to obtain simulations having approximately the same standard deviation of playing frequency. In practice we defined one α per note and dynamic (p, m f , f f ); its value ranges from 0.85 (D5) to 1.17 (B 3). During the optimization, the procedure defining ϕ for every bore relies on a maximal distance from the corresponding cloud of successful embouchures found during the preprocessing, above which a virtual embouchure is discarded. In the interest of keeping the length of this paper reasonable, we point the interested reader to the manual accompanying the source code repository for a detailed description of the procedure (framagit.org/rtournem/BrassOptimUsingSounds).

It is important to mention that for numerical reasons, the quantities F(x, ϕ i ) and SC(x, ϕ) are average values across a finite set of embouchures ϕ, randomly chosen and selected by the preprocessing. The consequence is that the objective function J(x) is non-deterministic, i.e., different objective function values may be obtained for the same x. In practice, a set of 1000 embouchures are simulated per note per dynamic range in order to keep the standard deviation on the numerical Brass Instrument optimization estimation of J(x) as low as possible according to the law of large numbers. This choice of 1000 embouchures is validated a posteriori given the small error bars in Figure 5. Additional "blackbox"

properties of the objective functions under consideration include:

• The evaluation of J(x) may fail due to difficulties to simulate notes (find virtual embouchures).

• It is not possible to reliably compute the gradient of J(x) because of the random selection process in the selection of the virtual embouchure.

• The evaluation of the objective function can be computationally expensive (between 3 and 20 minutes depending on the processor).

• We cannot assume smoothness of the objective (or constraint) functions.

To address these issues, we resort to the use of derivative-free optimization algorithms and a surrogate-assisted modeling strategy.

IV. SURROGATE-ASSISTED DERIVATIVE-FREE OPTIMIZATION

We use the rigorous derivative-free mesh adaptive direct search (MADS) optimization algorithm, which has convergence properties [START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF] and has been implemented in the NOMAD software package [START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm[END_REF]. Every iteration of the MADS algorithm consists of two steps: the optional search and the mandatory poll. The search step can implement any user-defined strategy to obtain promising candidates. The poll step determines candidates around the incumbent solution; it ensures the convergence of the algorithm towards a local optimum. Our strategy in the search step is to formulate and solve a surrogate problem to obtain a promising candidate, i.e., we use surrogate models of the computationally-intense simulation procedure to evaluate the objective and constraint function values. We then evaluate the real potential of this promising candidate using the physics-based simulations. In addition, when the MADS algorithm needs to proceed to the poll step, we use the surrogate models to rank-order the poll-generated candidates and then evaluate them opportunistically using the physics-based simulations. In this manner, we generate a large amount of information using computationally inexpensive surrogate models but make algorithmic decisions using the high-fidelity simulations.

More details are provided in the next sections.

A. Mesh Adaptive Direct Search

At each iteration k of the MADS algorithm, the trial points must lie on a mesh M k . The mesh size ∆ m k depends on the iteration number k and gets smaller as the optimization converges.

During each search step, a surrogate model Ĵ is built using previous evaluations of the objective function J. Then, a second instance of MADS is used to obtain the design that minimizes Ĵ.

This candidate design is then projected on the mesh M k and J is evaluated. If this candidate leads to an improvement of the solution, the surrogate model Ĵ is updated and the search is repeated.

Otherwise, the algorithm continues with the poll step. Two possible surrogate modeling techniques are described in the next section.

During each poll step, a set of candidates P k is generated on the mesh M k . The distance between the incumbent solution and the candidates P k is controlled by the poll parameter ∆ p k which, as ∆ m k , gets smaller as the optimization converges. The interested reader can refer to [START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF] for details. As mentioned earlier, we first rank-order the points of the set P k using the

V. APPLICATION

We consider three problem with ten design optimization variables: Int, SC Int and SC Dyn.

We have also considered these problems with two design variables similarly to [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF], but we ommit them to keep the paper length reasonable. Since it can be useful to optimization novices, we have made it available online at this link.

The initial bore x 0 is close to an existing trumpet bore (whose internal diameter has been measured with different balls of decreasing diameter inserted in the trumpet, and a gauge to measure their position inside the instrument).

The design problems are solved with a budget of 200 function (or blackbox) evaluations. This maximum number of function evaluations is defined empirically by observing the evolution of the objective functions in order to minimize computational cost (i.e., avoid unnecessary evaluations that do not improve the function value significantly, see Figure 4). Moreover, to ensure a reliable quantification of the efficiency of the optimization method, each problem is solved 10 times (with different starting points) for each of the 2 surrogate modeling approaches (using either an ensemble of surrogates according to [START_REF] Audet | Order-based error for managing ensembles of surrogates in derivative-free optimization[END_REF] or locally weighted regression models according to [START_REF] Talgorn | Locally weighted regression models for surrogate-assisted design optimization[END_REF]). Given that 3 design problems are studied, there is a total of 60 optimization jobs (3 problems × 2 surrogate approaches × 10 starting points).

In order to undertake this considerable computational endeavor, we relied on a high performance computing cluster, parallelizing the sound simulations on 2 Haswell Intel R Xeon R E5-2680 v3 2.5 GHz processors providing 24 cores for each of the 60 jobs. One evaluation of the objective function took between 2min 40sec and 3min 50sec for Int and SC Int (4000 sounds simulated), depending on the difficulty to simulate sounds for the considered geometry x. We provide sounds of the studied notes of the initial trumpet as supplementary materials 1 . These sounds are made of the first 6 harmonics of the permanent regime (no transients).

A. Design Optimization Problem with 10 Variables (10d problem)

The ten design variables represent the geometry of the leadpipe which is an important part of the bore that connects the mouthpiece to the tuning slide. The leadpipe, roughly conical, has a significant influence on the intonation and timbre of the instrument [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF].

Eleven parts of equal length (l=20mm) are considered. The design variables are the inner radii of the leadpipe at the connection between two parts (10 variables out of 12 control points because the initial and last control points are fixed at 4.64 and 5.83 mm, respectively). These 10 inner radii values span from 4.5 to 6 mm. The rest of the instrument corresponds approximately to a standard trumpet.

In this realistic design space, the high dimensionality requires efficient optimization: a discretization of the space with a granularity of 0.5 mm (20% of the range of each dimension) would necessitate 4 10 function evaluations (more than 1 million, which is not tractable in a reasonable computation time).

Figure 4 presents the performance of the optimization approaches for the three 10-d problem formulations. The initial value for J 1 is 8.7 cents. The initial value for J 2 and J 3 are 2.23 and 0.04 (value of SC), respectively, as can be seen in Figure 5.

On average, the LOWESS surrogate approach yields slightly better designs for J 1 (x, ϕ) and better results for J 3 (x, ϕ). For Int (J 1 ) the best objective function value obtained is 0.1 cents, For SC Dyn (J 3 ) the best objective function value obtained is 0.14 (value of SC), and the spectral centroid difference improvement is 0.1. This optimal bore improves 35 times the J 3 performance of the initial bore (0.14 / 0.04). Yet, measurements on an expert trumpet player playing B 4 concert pitch show a SC increase around 0.8 between a p and f f note for the 6 first harmonics, which is more than 5 times superior to the simulated values of the optimum. These results are in accordance with [START_REF] Petiot | Comparison of trumpets' sounds played by a musician or simulated by physical modelling[END_REF]. Even if the method optimizes the spectral centroid variation, the model does not allow the prediction of realistic spectral centroid values. Non-linear propagation in the bore should be taken into account.

Figure 5 provides a finer acoustical analysis of these results showing the contributions of each note to the objective functions for the initial bore and the optimal bore. The error bars are computed following the guidelines of the Joint Committee for Guides in Metrology [START_REF] Jcgm | Evaluation of measurement data -guide to the expression of uncertainty in measurement[END_REF]. For Int, the B 3 and F4 show significant improvements of 7.5 and 13.6 cent, respectively, which are supe- rior to the classical just-noticeable difference (JND) of 5 cents. B 4 being the relative reference, it is always considered perfectly in tune. Since the initial value of D5 is already low (4.8 cent) the improvement of 4.7 cent is slightly lower than the JND. The sound simulation results concerning intonation are in agreement with the impedance peaks frequency values except for D5, for which the optimal trumpet seems less in tune when considering the impedance peaks frequency values.

For SC -Int, the highest improvement is 0.08 (value of SC) for B 4 which is at the same level than the JND of 0.1 reported in [START_REF] Jeong | The dependence of timbre perception on the acoustics of the listening environment[END_REF]. Consequently, for B 3, F4 and D5, the improvements seem negligible, rising to 0.04, 0.01 (error bars level), and 0.5, respectively. For SC Dyn, the improvements for B 3 and F4 are above the JND (0.13 and 0.15, respectively), contrary to B 4 (0.01).

Figure 6 shows the leadpipes yielding the best value for each design problem.

As in [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF], the optima are counter intuitive since the leadpipe does not have a positive slope along the whole trumpet axis. This kind of shape for a leadpipe is not common among trumpets because it is very difficult to manufacture. The optimization algorithm was able to explore the design space in order to find unusual designs. Finally, the geometrical differences between the 3 optima and the initial leadpipe are on the order of the millimeter for several control radii, which will lead to noticeably different instruments when manufactured.

B. Discussion

Concerning intonation, the objective function value based on sound simulations is always below the one obtained only with the input impedance peaks (cf. black dots of Figure 5). This has been verified on more examples [START_REF] Tournemenne | Optimisation d'un instrument de musique de type cuivre basée sur des simulations sonores par modèle physique[END_REF], and would mean that the musician always plays

Brass Instrument optimization

Regarding the surrogate modeling approach, it is hard to tell whether one is consistently better than the other as they seem to exhibit the same exploration to exploitation ratio of the design space. It should be noted that the locally weighted regression models take longer to compute than the ensemble of surrogates in high-dimensional problems (in our case, the difference is one hour vs. few minutes for ten variables).

VI. CONCLUSIONS

In this paper, we extended a new paradigm for design optimization of brass instruments. This new paradigm goes beyond impedance based optimization. (Braden et al., 2009; Guilloteau, 2015; Kausel, 2001; Macaluso and Dalmont, 2011; Noreland et al., 2013 , 2010; Poirson et al., 2007).

We consider the optimization of objective functions (possibly subject to constraints) based directly on the sounds spectrum, which, to the best of our knowledge, is a novel approach. The originality of the approach lies in the fact that the objective function is not limited to a characterization of the instrument alone, but includes virtual musicians in an physical model, to optimize directly the instrument sounds. The main contribution of this paper is to demonstrate how physics-based sound simulations can be integrated in an iterative optimization algorithm, which requires that simulations converge automatically toward auto-oscillations for every considered bore of the design space, without any assistance of the user. A second contribution concerns the optimization of objective functions based on a particular dimension of the timbre of the sounds, the spectral centroid. Applied to the optimization of a trumpet, the results show that the optimization method, based on the MADS algorithm, is efficient to define optimal solutions in a reasonable computation time, with or without constraints, for problems up to 10 design variables.

While the approach shows promising potential, there is room for improvement. Even if it is efficient to reproduce differences between instruments concerning intonation and spectral centroid, the elementary model could be improved to generate a more realistic sound spectrum. The sound optimized in this work is the sound in the mouthpiece. Even if this does not change the principle of the method presented, it could be interesting to define a relevant radiated pressure outside the instrument, and to include it in the optimization considerations. Another valuable contribution to this numerical study concerns the manufacturing of the optimal instruments and their objective and subjective study. This would help evaluate actual improvement. Regarding implementation, the influence of the selection process of the virtual embouchure could be investigated further, as it may provide more robust descriptors of the ease of playing of the considered instruments. Regarding the methodology, a study of temporal sound simulations could lead to new classes of objective functions, such as attack times. More ambitious still, the inclusion of non-linear propagation in simulations would produce more realistic sounds improving the objective functions considered in this work. Finally, the main challenge in the design of musical instrument lies in the definition of judicious objective functions providing actual insight of the instrument intrinsic quality.

This task may be accomplished by working side-by-side with instrument makers whose empirical understanding may be translated to computational principles and models.

FIG. 1 .

 1 FIG. 1. Simulated input impedance Z of a B trumpet (magnitude), highlighting resonances 2, 3, 4, and 5

  FIG. 2. Representation of the outward striking model of the lips, with the definition of the variables of the

Figure 3

 3 Figure3describes the flowchart of the process for optimizing the shape of a trumpet bore using

  FIG.6. Representation of the leadpipe inner radius along the instrument axis; the black dotted line to the

TABLE I .

 I Values of the control parameters for the simulations considered in the study (virtual musician

	embouchure)
	Parameter	Symbol (units) Value
	Resonance frequency of the lips	f (Hz)	130 to 480
	Mass per area of the lips	µ (kg/m 2 )	1 to 6
	Pressure in the mouth	P m (kPa)	1 to 12
	Width of the lips	b (mm)	10
	Rest value of the opening height	h 0 (mm)	0.1
	Quality factor of the resonance	Q	3

  On the top part, the columns represent from left to right J 1 value and the non zero ETDs. Two bores are evaluated, the initial bore of all the optimization jobs and the best bore from the most successful job (out of the 20 available). The error bars equal to 2 standard deviations of the estimated quantity. Equivalent

							SC -Int
	i (in cent) J 1 or ETD	-30 -20 -10 0 10 20	8.7 17.2	0.1 10	-7.7 -26.4	0.2 -8.3	-13.7 -21.7	-0.1 -12.9	4.8 3.5 initial 0.1 8.9 optimal impedance
				J 1		ETD 2	ETD 3		ETD 5
							SC -Int
	2 or SC i (SC value)	2.1 2.3 2.5	2.23	2.27	2.45	2.49	2.33	2.34	2.14	2.22	1.99 initial 2.04 optimal
	J	1.9
			J 2	SC 2	SC 3	SC 4	SC 5
							SC Dyn
	J 3 or ∆SC i (SC value)	0 0.05 0.1 0.15 0.2		0.04	0.14	0.01	0.14	0.04	0.19		0.07 optimal 0.08 initial
				J 3	∆SC 2	∆SC 3		∆SC 4
	FIG. 5.

plots for the design problems SC -Int and SC Dyn are drawn below. SC i is SC restricted to ϕ i . In addition, the values of J 1 and the corresponding ETDs, obtained with the input impedance peaks frequency values, are presented in black dots.

ACKNOWLEDGMENTS

The authors are grateful to the NOMAD team (Prof. Charles Audet, Prof. Sébastien Le Digabel and Dr. Christophe Tribes), Prof. Saïd Moussaoui for his help with the optimization problem formulation, Prof. Jean-Pierre Dalmont for his precious help with the input impedance models,

surrogate model Ĵ. The physics-based sound simulation model J is then used to evaluate the points of the set P k using an opportunistic strategy: If a point is feasible and leads to an improvement of the objective function, the evaluation process is aborted and the algorithm iterates. Note that if a more feasible point is found during the poll step, the mesh and poll parameters are increased so that the algorithm can explore other areas of the design space. Otherwise, these parameters are reduced, which means that the iteration will operate in a closer neighborhood of the design space.

B. Surrogate Modeling Strategy

We consider two surrogate model approaches in this study:

• An ensemble of surrogate models approach where several surrogates of different types and with different modeling parameters are built and updated while selecting the model that fits the data best at each iteration. We use the same ensemble of surrogates as in [START_REF] Audet | Order-based error for managing ensembles of surrogates in derivative-free optimization[END_REF] and [START_REF] Tournemenne | Brass instruments design using physics-based sound simulation models and surrogate-assisted derivativefree optimization[END_REF]: 6 polynomial regression models, 5 kernel smoothing models, and 6 radial basis function (RBF) models. The best model is chosen according to the order-error with cross-validation (OECV) metric presented in [START_REF] Audet | Order-based error for managing ensembles of surrogates in derivative-free optimization[END_REF].

• A LOcally WEighted Scatterplot Smoothing [START_REF] Talgorn | Locally weighted regression models for surrogate-assisted design optimization[END_REF] (LOWESS) surrogate model. This model consists of building a local polynomial regression around the point x where we wish to predict the objective function. In the construction of this local regression, data points close to x are given more importance than those further away from x.

Society of America 97(3), 1929-1939. 573


Marion Maisonobe
email: marion.maisonobe@cnrs.fr

Laurent Jégou
email: laurent.jegou@univ-tlse2.fr

Nikita Yakimovich

Guillaume Cabanac
email: guillaume.cabanac@univ-tlse3.fr

NETSCITY: a geospatial application to analyse and map world scale

published or not. The documents may come

Introduction

Practitioners in the field of bibliometrics contributed a variety of software tools to simplify the processing and mapping of bibliographic data. We can think of the online dashboards available through the Web of Science Core Collection and Scopus websites. There is also a variety of free bibliometric mapping software tools such as VOSViewer (van Eck & Waltman, 2010) and CiteSpace [START_REF] Chen | CiteSpace II: Detecting and visualizing emerging trends and transient patterns in scientific literature[END_REF]. According to the review proposed by [START_REF] Cobo | Science mapping software tools: Review, Analysis, and Cooperative Study Among Tools[END_REF], which focuses on these latter tools together with Bibexcel, CoPalRed, IN-SPIRE, Leydesdorff's programs, Network Workbench Tool, Sci 2 Tool, and Vantage Point, existing software tools offer four different types of bibliometric analyses: burst detection, geospatial analysis, network analysis, and temporal analysis. Among the nine tools reviewed by Cobo et al., network analysis is the most prominent one. Indeed, existing tools primarily allow computing and/or mapping of bibliometric networks of words, authors, and journals derived from co-occurrence data. Geospatial analysis, on the contrary, is the less common type of available analyses. It is included in only three of the reviewed software tools: CiteSpace, Sci 2 Tool, and VantagePoint. The first two allow network data visualisation over a world map using Google Earth Maps or Yahoo! Maps. They also include geocoding capabilities but the choice of the spatial resolution is only dual: either the data are clustered at the country level, or they are clustered at the street levelstreet addresses used in the publications are then converted to coordinates in latitude and longitude [START_REF] Chen | CiteSpace: A Practical Guide for Mapping Scientific Literature[END_REF]. In addition to geocoding and mapping capabilities, NETSCITY, the online application for geographic analysis of research output that we present in this article, features an intermediary level of geographic aggregation, which is the urban area level. Actually, NETSCITY tackles four main issues in the business of mapping bibliographic data: the geocoding accuracy, the geographic aggregation, the counting/fractioning issue, and the mapping issue. NETSCITY is a free geoweb application developed to promote the spatial scientometrics method designed by our team to research on the world geography of scientific production and collaboration. Over the past 10 years, our team has been working on geocoding, clustering, and mapping the contemporary geography of scientific activities using bibliometric data. To measure production share and collaboration intensity, we delineated urban area perimeters covering geocoded addresses. This aggregation step proved necessary to work with statistically comparable geographic entities at the world level, which administrative municipalities or street addresses are not. While sharing a set of 495 urban area perimeters in an open access publication (Maisonobe et al., 2018) we have developed NETSCITY, which allows any user:

-to upload a set of references to scientific materials; -to geocode them; -to aggregate them at the urban area and the country levels; -to compute the number of publications per spatial unit and the number of collaborations between spatial units using various counting methods; -to visualise these variables both on a geographic and on a network map. At each step, the user is free to download the results of the current process. Some users can be interested in the results of the geocoding process only, others in the results of the aggregation process, others in the results of the counting process, and others in the results of the mapping process. NETSCITY thus fulfils a various number of spatial analytics purposes. We designed it is as a key application to ensure a better reproducibility in the spatial bibliometric field [START_REF] Frenken | Spatial scientometrics: Towards a cumulative research program[END_REF][START_REF] Cobo | Quantifying the reproducibility of scientometric analyses: a case study[END_REF]. We also plan to release the source code of the application.

In this paper, we first discuss the inputs of NETSCITY for the spatial bibliometric field. Then, we describe the inner workings of NETSCITY: from upstream geocoding down to mapping processes. Next, we show the results of NETSCITY when applied to a specific set of references. We conclude with future works involving additional levels of geographic aggregation, counting methods, and mapping options.

Mapping the geography of science

Geocoding

Mapping scientific activities is a five-century-old craft, with seminal maps of science showcased in Börner's (2010) Atlas of Science. In 2010, Loet Leydesdorff and Olle Persson (2010) published a landmark article discussing the opportunities of geocoding and mapping Web of Science and Scopus data at the city and institution levels. One year before, [START_REF] Frenken | Spatial scientometrics: Towards a cumulative research program[END_REF], were identifying a spatial turn in scientometrics that has been confirmed later [START_REF] Frenken | Spatial Scientometrics and Scholarly Impact: A Review of Recent Studies, Tools, and Methods[END_REF]. The improving capabilities and availability of online geocoding tools, the trend toward territorialised policies of science, and the renewed enthusiasm for global benchmarking of cities and higher education institutions prompted this spatial turn. At the same time, our research team, mainly based in Toulouse (southwestern France), identified a need for strengthening the ties between bibliometrics and geomatics. Processing the spatial information included in bibliometric datasets by relying only on the results of online geocoding tools proved insufficient, such as Google Maps API [START_REF] Jégou | Toward spatially referenced academic data at global scale: the full geocoding of Wos-Datasets, methods and results[END_REF]. Many errors result from trying to geocode the entire address strings included in bibliometric datasets. Let us consider the following example: an author of a 2009 publication reported the following address: "Chang Gung Mem Hosp, Tao Yuan, Taiwan". It corresponds to the Chang Gung Hospital, in Taoyuan District, Taoyuan county of Taiwan Island, 15 km west of the capital, Taipei. Google Maps API hesitates between the hospital itself and the Chang Gung Memorial metro station. Other web geocoding services confuse this address with several locations in the district. Indeed, the district name is similar to the county name and to the island name.

If the required level of aggregation is the urban level, we proved that it is more relevant to begin by spotting the portion of the address string corresponding to the city, the province / state and the country. The latter is the easiest to find since it always ends the string. As shown in the aforementioned example, tagging the city and the province is a more complicated task because their order of appearance might differ from one address to another and the province is not always specified (non-federal countries such as France do not use them in postal addresses). Even by pre-processing the data to ensure better geocoding results, some errors or missing values still occur: some of them result from homonyms, which require additional spatial information, other result from misspelled toponyms. Here we should highlight that geocoding as we performed it is a semi-automatic process: a geographer expert reviews the results, focussing on the records affiliated to the most prominent scientific cities or the most prone to potential errors (we devoted additional work on addresses having undergone an alphabet transliteration and/or with postal system specificities like Indonesia or Taiwan). After years of geocoding work on the entire contents of the Web of Science Core Collection and on Scopus datasets, our research team improved the geocoding results one obtains when using these databases [START_REF] Jégou | Toward spatially referenced academic data at global scale: the full geocoding of Wos-Datasets, methods and results[END_REF]. Capitalizing on this previous work, the geocoding tool implemented in NETSCITY detects and corrects much of the common misspellings found in Web of Science and Scopus data (by simplifying the given addresses, comparing them with a list of known misspellings variants and, in the absence of a match, by using internal and external online geocoding services) -these stem from authors, publishers, typists working for such database vendors. For the remaining errors, NETSCITY users can amend the records as they see fit (see next section).

Geographic levels of analysis

To study the world geography of science, our team strived to determine an adequate urban level of data aggregation. For studying the geography of science before 2010, most scholars were relying on the country level and only a few studies tackled it at the urban level. The latter has attracted a growing interest since the 2000s following city and regions' empowerment and the multiplication of city rankings [START_REF] Bornmann | Spatial bibliometrics on the city level[END_REF]. Two types of studies can be distinguished: those focusing on a limited number of urban regions (e.g., [START_REF] Matthiessen | World Cities of Scientific Knowledge: Systems, Networks and Potential Dynamics. An Analysis Based on Bibliometric Indicators[END_REF][START_REF] Hoekman | Research collaboration at a distance: Changing spatial patterns of scientific collaboration within Europe[END_REF][START_REF] Nomaler | On Scaling of Scientific Knowledge Production in U.S. Metropolitan Areas[END_REF]; and those encompassing the publishing localities of the entire world (e.g., [START_REF] Waltman | Globalisation of science in kilometres[END_REF][START_REF] Pan | World citation and collaboration networks: uncovering the role of geography in science[END_REF][START_REF] Csomós | A spatial scientometric analysis of the publication output of cities worldwide[END_REF]. In the first case, the authors tend to reuse existing sets of administrative perimeters: United States MSA, European NUTS and, in the second case, the authors do not tend to aggregate the geocoding localities. To get a global panorama of the world scientific production, neither of these two approaches is satisfying. Existing sets of urban area perimeters fulfil specific purposes (whether they are defined for a national or a continental scope, or they are limited to the most populated areas only). Combining existing sets in order to encompass the entire world production might be a solution but it leads to comparative biases since the definitions used to produce these various sets are diverse. Not aggregating the data at all also leads to comparative biases since municipalities' boundaries depend on very heterogeneous criteria differing from one country to another (Maisonobe et al., 2018). For example, not aggregating the data will lead to count apart the numerous publications authored from Bethesda (Maryland, USA) of those of Washington City (DC, USA). Bethesda is a suburb of Washington, which include important medical institutions such as the NIH. To overcome these issues our team designed a unique set of urban area perimeters encompassing all publishing localities identified after geocoding the entire 1999-2014 contents of the Web of Science Core Collection. Our urban delineations consider the distribution of the world population density and the Euclidean distance between publishing localities. So far, NETSCITY allows its users to aggregate their data at the level of this set of urban area perimeters as well as at the country level.

Counting method

The next issue addressed by NETSCITY is the counting conundrum. Since about one third of the entire world publications are authored from more than one urban area, one needs to decide how to measure urban areas' contribution to science production. There exist many ways of assigning a number of publications to statistical units. [START_REF] Gauffriau | Comparisons of results of publication counting using different methods[END_REF] gives a rich overview of existing methods and highlight the need for more transparency on the methods used since the choice of a counting method significantly influences research findings in bibliometrics. When focusing on the geography of science, the elementary unit of analysis can be the author or the institution, but given the biases of existing sources, it is more often the address. Indeed, in pre-2008 bibliometric records, the list of addresses is not always linked to the list of authors. As a result, it is possible to derive a number of authors per publication as well as a number of addresses but not a number of authors per address per publication. In addition, authors' addresses do not necessarily correspond to authors' institutions since an address may refer to several institutions. Conversely, different addresses may refer to the same institution with varying research teams or university departments. Unless pre-processing the data accordingly, the postal address is the most elementary counting unit of a bibliometric dataset for geographic purpose. As a result, there are two main ways of measuring the scientific production of an urban area: adding up the number of addresses per urban area involved in the publication or considering the number of different urban areas involved in the publication. Similarly, to measure the scientific production of a country, we can consider three different values: the number of addresses per country involved in the publication, the number of urban areas per country involved in the publication, or the number of involved countries. In some disciplines such as in chemistry, only the first address or the corresponding address might be relevant. Other addresses would not be counted in this latter case. In addition to the choice of a counting unit, the counting issue requires to arbitrate between full and fractional countings [START_REF] Van Hooydonk | Fractional counting of multiauthored publications: Consequences for the impact of authors[END_REF]. In the first case, we consider the total number of addresses/urban areas/countries while in the second case, we fraction the credit of the publication (i.e., one unit) so that the sum of each fractioned credit total one. Since fractioning avoids counting a single country contribution multiple times, it is the most preferred method in bibliometrics. As shown by [START_REF] Leydesdorff | Full and Fractional Counting in Bibliometric Networks[END_REF], the choice of a counting method applies both to production data (computing a number of publications per spatial entity) and to collaboration data (computing a collaboration weight between co-authoring entities). With NETSCITY, it is possible to control for the computation of both indicators.

Mapping issue

Mapping issues addressed in bibliometric research are traditionally limited to the visualisation of co-occurrence networks in a relational space. Waltman & van Eck (2010) distinguish two main types of bibliometric maps, namely graph-based maps and distance-based maps. The former refer to mapping the presence or absence of a relation between entities while the latter refer to mapping the distance between entities according to the intensity of their relations. There exists a large range of similarity metrics to compute a relational distance between bibliometric entities: cosine, association strength, Jaccard index, Pearson index and, of course, the raw cooccurrence number (either full or fractional as shown in the previous sub-section) [START_REF] Boyack | Mapping the backbone of science[END_REF][START_REF] Van Eck | How to Normalize Co-occurrence Data? An Analysis of Some Well-Known Similarity Measures[END_REF]. When used to produce a distance-based map, similarity metrics can be selected according to different mapping algorithms. Kamada Kawai and Fruchterman Reingold are the most commonly used techniques in the field of bibliometric as well as in the broader field of network analysis. To adapt to power-law characteristics of most bibliometric distributions, [START_REF] Waltman | A unified approach to mapping and clustering of bibliometric networks[END_REF] propose a unified way of mapping and clustering bibliometric networks. VoSViewer as well as other network analysis software such as Pajek include this method. Compared to this latter issue, mapping issues related to the production of geographic maps are mostly overlooked. The world maps used in spatial bibliometrics are often screenshots taken from Google Earth or Yahoo! Maps' base maps. Instead of benefitting from the capabilities of thematic cartography to map production and collaboration data, there seems to be a preference for online and interactive visualisation. However, we believe the two types of graphic representation -static and interactive -should be valued and different visualisation methods must be chosen accordingly, within the context of geomatics. In this respect, the literature in GIS and cartography about flow maps and the development of geoweb applications might be of relevance [START_REF] Dodge | Geographic Visualization: concepts, Tools and Applications[END_REF]. While NETSCITY displays interactive geographic and network maps, it also allows the downloading of all underlying data necessary to generate static maps. This feature of the application is about to include more adjustable parameters, interactive linked graphs, and a vector graphics output. In forthcoming developments, NETSCITY will also give access to the Netmap visualisation interface1 allowing the simultaneous exploration of the geographic and the network maps: by clicking on an item on the graph, the item and its relations are highlighted on the map and vice versa (Maisonobe & Jégou, 2018). As demonstrated in [START_REF] Andurand | Les mondes savants et leur visualisation, de l'Antiquité à aujourd'hui[END_REF] and [START_REF] Bach | NetworkCube: bringing dynamic network visualizations to domain scientists[END_REF], combining several types of graphic representation enriches the data exploration and analysis experience.

In what follows, we explain how to use the NETSCITY geoweb application and we showcase results obtained by applying NETSCITY to a set of references extracted from the Web of Science.

Features implemented in NETSCITY

In this section, we provide a description of the geoweb application from the data import and processing step to the visualization step. We also discuss technical implementations and directions for improvements throughout the section. NETSCITY is designed for different information needs and types of users:

-The researchers interested in mapping their own field of research or their own publication record. -The librarians interested in organizing sets of references and mapping the geographic relations between scientific references. -The policy officers or international experts interested in scientific evaluation.

-The students or journalists interested in mapping a field of knowledge or the production of a specific institution, author, or place. -The science studies researchers interested -as we are -in mapping the geography of science.

The home page of NETCITY, presenting the project, is available following this web link: https://irit.fr/netscity The main hub features four types of actions are available: 1) import a set of references extracted from the Web of Science or Scopus online websites, or personal CSV records, 2) explore the results (lists or data-visualizations) 3) export the results 4) visit the help pages.

Importing a dataset from Web of Science, Scopus or else

In the current version, users can only import one type of Web of Science file. From the Web of Science web interface, they must download a set of references using the following method: 1) Search for a list of records to export and select "Save to Other File Format", 2) Export all records returned (up to 500) by entering "1" to "[the total number of records]", 3) Select "Full Record and Cited References" to include addresses and times cited information, 4) Select "Tab-Delimitated (Mac/Win, UTF-8)" from the "save" drop-down menu. Similarly, for a Scopus file, the procedure is the following: 1) Click on the Export icon2 , 2) From the pop-up menu, select "CSV" format for your file format and "Citation" as well as "Bibliographical information" as the information to export, 3) Click "Export", select "CSV -Only the first 2,000 documents" and click "Export" again.

It is also possible to import a personal CSV file. To do so, users must specify the name of the ID column and the Year column. The name of the columns storing the spatial information is also needed. Several options are possible since the spatial information can be stored:

-In entire postal addresses' strings, -In three pre-processed columns corresponding to the triples: "city", "province", "country", -Both in entire addresses' strings (including the institutions) and in pre-processed columns.

The text file can be compressed before the upload, to speed up the transfer. Once one file is uploaded in NETSCITY, it is possible to import another file without deleting the previously uploaded one. This allows us to deal with a bigger dataset and adding the data progressively, coming from various bibliometric sources. According to users' needs, we are eager to improve the import options and open it to other types of sources and formats.

Geocoding and using the error-correction tool After the import step, the user gets a report detailing the quality of the geocoding process. In Figure 1, the geocoding process has been applied to a homemade CSV file. The bibliometric records included in this file are coming from different sources: Web of Science, Scopus and a French open archive (HAL). As a result, the quality of the geocoding process is poorer than when applied on Web of Science or Scopus files only. However, nearly all the 1,747 addresses included in this file were geocoded. The majority were geocoded with no recourse to an API 3 , which means that our internal geocoding knowledge base was sufficient to resolve the geographic coordinates. The remaining addresses were geocoded using either a state-of-the-art gazetteer (Geonames), integrated in the application, or an online geocoding API (LocationIQ service).

Figure 2. The error-correction interface allowing users to amend the results of the geocoding process

As for non-geocoded addresses, NETSCITY offers the possibility to type some missing spatial information or corrections to complete the geocoding process (Figure 2). Verifying or filling the content of the fields in red (i.e., "city", "province" and "country") enhances the geocoding process. If the geocoding fails despite these complements or corrections, it is still possible to add manually the geographical coordinates corresponding to the non-geocoded addresses. In a future version of NETSCITY, we plan to let users modify the result of the geocoding process for geocoded addresses to thwart false positives. After this first step, users can export the result of the geocoding process, add another file or directly go to the data exploration tools. When choosing to export the results of the geocoding process, NETSCITY generates a list of addresses together with the eventually corrected name of the corresponding cities, provinces and countries (these are sometimes misspelled), the publications ID and the geographic coordinates in a CSV file.

Exploring the dataset through statistical tables and maps Four types of data exploration are available. The users are presented with:

-Production data through a map and a table listing the number of publications per spatial entity (urban area/country) -Collaboration data through a map and a table listing the number of collaborations between all pairs of spatial entities (urban areas/countries) When accessing these different views, the users can change the counting method and unit so that they can immediately see the effects of their choice on the ranking (table view) and on the visualisation (map view).

In the current version, two levels of counting units are available: the users can opt for 1) counting the number of addresses per spatial entity (urban area/country) or 2) counting the number of urban areas or countries. Then, the users are free to choose between full and fractional counting methods. Applied to collaboration data at the urban area level, the fractional counting method implies that by adding the weights of all the links, we obtain the total number of co-written publications of the dataset. It implies weighting the links according to the number of interurban links per co-publication. For instance, if a given publication stems from three different urban areas, each inter-urban link receives 1/3 as a weight for this publication. More generally, if a publication is co-signed from 𝑛 urban areas, each pair of urban areas (A, B) with A < B is assigned a value 𝑙 equals to: 1/𝑛(𝑛 -1)/2 = 2/(𝑛(𝑛 -1)). When exploring the numerical data online, all temporal information is overlooked. Longitudinal analysis is available at the export stage, however: the option "with the evolution by years" computes the different variables (number of publications and number of collaborations) on an annual basis. Then, the export final step allows, on the one hand, opting for full or fractional counting and, on the other hand, selecting between CSV or JSON formats.

Availability of NETSCITY as an open source application

The source code is to be released in an online open-source forge (like GitHub) with a suitable license (like one compatible with the well-used GNU-GPL), after several enhancements and another pass of testing and debugging. The project is still ongoing, and we are eager to collect feedback and ideas to expand its capabilities, to tune it to new needs, and to enhance interoperability with other analysis applications such as the CorTexT project4 or VOSViewer. Together with several bibliometrics tools, the geocoding tool and the aggregation tool included in NETSCITY could be available as off-the-shelf modules.

Geographical underpinning of a set of bibliographic references: a case study using NETSCITY

To test NETSCITY on a specific set of publications, we decided to focus on the publications referring in their titles, abstracts or keywords to the use of a common device: ROV or AUV, which is a submarine robot. Therefore, we issued the following topic search in the Web of Science Core Collection: TS= ("ROV" OR "AUV") AND "robot*" AND "underwater", all publication years combined. This query returned 1,004 records published between 1991 and 2017 (Figure 3). We consider this WOS result in the running example that follows.

Figure 3. Result page for a sample Web of Science query with 1,004 records

Although the Web of Science interface cannot reveal the leading urban areas producing knowledge on the matter under consideration, NETSCITY lists the top 10 urban areas which publish the most about ROV/AUV (Table 1) at once. Then, we can identify the top 10 most intense interurban collaborations about ROV/AUV according to a different counting method (Table 2). In addition, we can interactively explore the world production map and the world collaboration map per urban area and per country using different counting methods. Figure 4 and Figure 5 are screenshots taken from NETSCITY cartographic zoomable views. From these four views (we could also have chosen to discuss the networks views available in NETSCITY), we can notably observe that the ROV/AUV topic involves publications coming from all over the world, but mainly from maritime cities (Figure 5). We can also notice an important Italian cluster of scientific collaborations on this research issue (Table 2). It is worth noting here that we obtained all these results with a few clicks only. NETSCITY constitutes an effective application to efficiently and visualise and explore the geography of a bibliometric dataset according to different levels of analysis and counting methods.

Conclusion

The field of bibliometrics already develops and provides a variety of software to process bibliometric data without requiring any programming skills. However, most existing software fail to account for the geographical complexity of bibliometric data, which leads to misleading results. In addition, many applications generate graph-based maps of locations but, to the best of our knowledge, there is no fully integrated solution to parse and clean the affiliations recorded in a set of references, geolocalising them, clustering them at the requested level of analysis, and mapping them both on a world base map and in a relational space. The goal of NETSCITY is to offer and promote such a solution. Concisely, NETSCITY is a new application designed for the general public to perform spatial bibliometrics at the click of a mouse. It computes metrics and produces maps abiding by the state-of-the-art methods of geography of science. We wish that NETSCITY contributes to overcome the pitfalls of spatial data and better inform policy makers regarding territorialised policies worldwide.

Figure 1 .

 1 Figure 1. The result of the geocoding process showing a breakdown of the methods used to geocode the addresses from the uploaded addresses

Figure 4 .

 4 Figure 4. Production map for the 1,004 WOS records of the running example, after geocoding and aggregation are performed. Fractional counts at the urban area level.

Figure 5 .

 5 Figure 5. Collaboration map for the 1,004 WOS records of the running example, after geocoding and aggregation are performed. Fractional counts at the country level.

Table 1 . Top ten urban areas which publish the most in the ROV/AUV field (Web of Science)

 1

	Urban area	Country	Full number
			of
			publications
	Tokyo	Japan	44
	Boston	United-States	41
	Genoa	Italy	41
	Shanghai	China	39
	Girona	Spain	37
	Pisa	Italy	30
	Florence	Italy	28
	Woods-Hole	United-States	26
	Edinburgh	United-Kingdom	24
	Singapore	Singapore	23

Table 2 . Top ten most intense interurban collaboration in the ROV/AUV field (Web of Science)

 2

	Urban area 1	Country 1	Urban area 2	Country 2	Fractional
					number of
					collaborations
	Florence	Italy	Pisa	Italy	30.0
	Waterloo-Guelph	Canada	Shanghai	China	16.0
	Genoa	Italy	Pisa	Italy	14.0
	Florence	Italy	Genoa	Italy	9.5
	Geelong	Australia	San-Antonio United-States	9.0
	Oslo	Norway	Trondheim	Norway	8.5
	Kitakyushu	Japan	Tokyo	Japan	8.0
	Brest	France	Shanghai	China	7.5
	Lisbon	Portugal	Porto	Portugal	7.0
	Woods-Hole	United-States	Baltimore	United-States	6.5

Source code of the prototype: https://gitlab.com/ljegou/netmap Demo: http://www.geotests.net/test/net_map/science_mep.html

Depending on how you use Scopus, this might be displayed according to your export preference. To select a different export method than what is displayed, click on the drop-down arrow

An external service providing geocoding information.

https://www.cortext.net

Acknowledgments

We would like to thank our entire geography of science research team, and especially Béatrice Milard, Michel Grossetti, and Denis Eckert. This research is supported by LABEX SMS (ANR-11-LABX-0066) under project codenamed Netscience.