G Bourgeois

V Meli

R Antonelli

C Socquet-Clerc

T Magis

F Laulagnet

B Hemard

M Bernard

L Fellouh

P Dezest

J Krawczyk

S Dominguez

F Baudin

J Garrione

C Pellissier

J.-A Dallery

N Castellani

M.-C Cyrille

C Charpin

F Andrieu

G Navarro

Crossbar Arrays based on "Wall" Phase-Change Memory (PCM) and Ovonic-Threshold Switching (OTS) Selector: a Device Integration Challenge Towards New Computing Paradigms in Embedded Applications

Keywords: Crossbar, "Wall" PCM, OTS selector)

HAL is

Introduction

Phase-Change Memory (PCM) and Ovonic-Threshold Switching (OTS) selector have been successfully co-integrated in 3D stacked Crossbar arrays definitely reaching a high maturity, until the mass production of devices targeting stand-alone Storage Class Memory applications [1]. In order to enable neuromorphic computing architecture design, for example to speed up matrix-vector multiplication operations [2], Crossbar Resistive Memory is considered a valuable solution to be implemented in future embedded architectures, in order to support recognition and prediction tasks required in autonomous systems (e.g. driving). In this work, we demonstrate the co-integration of a PCM resistive device (1R) based on "self-aligned" (SA) Wall structure [3] with an OTS selector device (1S) by adding a "Cross" patterning (i.e. "Double-Patterned Self-Aligned" or DPSA integration) perpendicular to the first SA one. We address the issues inherent to our DPSA integration and in particular related to the bitline realization, which should reliably recover the contact over the 1S1R devices, and to the lithography strict alignment to allow device scaling down to tens of nm.

We show a good control of the effects of etching/stripping chemistry on the integrated materials by TEM/EDX analyses and finally we provide electrical results demonstrating successful programming operations in devices with a surface size down to 60 nm×80 nm.

Results and Discussions

A. PCM and OTS alloys

Our 1S1R DPSA structure is schematically represented in Fig. 1a. The PCM layer is integrated over a heater element (H), and it is followed by the OTS layer based on a GeSbSeN alloy [4], [5] (Fig. 1b). The intermixing between the two layers is hindered by an intermediate TiN electrode. We fabricated 1S1R devices integrating two different GeSbTe PCM alloys A and B, respectively with low and high Ge content. They were coupled with two OTS thicknesses (i.e. t1 and t2 with t1 < t2), leading to a total of four different PCM-OTS stacks.

B. DPSA 1S1R cell fabrication

We fabricated our 1S1R devices and arrays in the BEOL of LETI Memory Advanced Demonstrator based on 130 nm CMOS technology. Fig. 2 summaries the key steps of the DPSA process integration described by Coventor SEMulator3D  modeling platform. A first Wall (W) SA patterning (x direction) of the PCM-OTS stack involving also the heater definition is followed by a Cross (C) perpendicular patterning (y direction). The recovery of the contact on the top electrode after Wall patterning represents the first challenge in such integration. We compared two approaches: "Dry Etch-back" (Fig. 3a) and "Full CMP" (Fig. 3b) (i.e. chemical-mechanical polishing). The first one demonstrated (in SEM analyses) a high dependency on the devices pitch, with a consequent high variability. The second one, based on the selectivity of SiO 2 polishing wrt a SiN layer used as stop layer demonstrated a good reproducibility. In order to obtain fast prototyping of scaled devices in x and y directions (targeting devices down to 50 nm), we took advantage of Variable Shape Beam VISTEC  SB3054 lithography tool allowing alignments control thanks to modelling and die-to-die correction. Being the C/H alignment crucial for the device functionality (i.e. to obtain an active region perfectly centered into the PCM layer), we designed dedicated structures with a large Wall size and for different Cross sizes, connecting many devices in parallel for the electrical evaluation of the C/H alignment at end of fabrication. Fig. 4a shows the results obtained on a wafer integrating only a PCM layer (w/o OTS). The resistance increases with the decrease of the Cross size, likely related to a gradual loss of reliability on the C/H alignment. The inset reports the typical overlay (OV) obtained in a whole lot for both W/H and C/W alignments with a mean OV+3 < 30 nm for both. This is affecting the alignment for a Cross size below 60 nm (a solution for a direct C/H alignment is under validation). On the contrary, we observe the expected almost linear relation between the devices virgin resistance (R V ) and the reciprocal of the Wall size (Fig. 4b). TEM analyses performed on a 60×60 device (i.e. Wall×Cross sizes are expressed in nm) integrating the PCM-OTS stack (Fig. 5a) highlight that alignment control is mandatory to target small sizes functionality. Indeed, we observe that H is not perfectly aligned with C in the device observed. EDX profiles (Fig. 5b) evidence a uniform distribution of the elements perpendicularly to the OTS/SiN interface. Indeed, chalcogenide materials are sensitive to halogens [6], then a step-by-step etching strategy with dedicated chemistry recipes for each step was applied during the Wall etching (i.e. starting from the top electrode: TiN / OTS / TiN / PCM / heater) and Cross etching (i.e. TiN / OTS / TiN / PCM) to prevent from OTS and PCM layers damage.

C. DPSA 1S1R electrical characterization

The firing operation (i.e. initialization step) was investigated in the four PCM-OTS stacks, showing a decreasing of the firing voltage (V FIRING ) when reducing the OTS thickness from t2 to t1, and when reducing the Ge content in the PCM alloy from B to A (Fig. 6a). As expected, V FIRING does not depend on the size of the device (i.e. electric-field-dominated phenomenon) [7], as shown in the results that are equivalent for 60×60 and 300×300 devices. Device functionality was statistically verified in 1 kb arrays of transistor-selected 1S1R devices (1T1S1R). As an example, in Fig. 6b we report the array map of the current after the devices activation, meaning the dynamic current measured during the pulse application. Only four devices are not switching (< 0.4%).

Threshold Voltage (V th ) versus programming current plots of Fig. 7 show the typical SET-RESET characteristics for 1S1R devices. We compared the behavior of 60×80 and 300×300 cells, both showing a reading window of about 1.5 V, compatibly with the V th of Ge-rich GeSbTe devices [8] (i.e. the 1S1R max voltage reading window is equal to the PCM V th [9]), with well distinguishable SET and RESET programming current regions. Current reduction in 60×80 cells is in line with the reduced heater surface wrt 300×300. SET/RESET programming operations were performed in 1T1S1R 1 kb arrays based on 60×80 devices along 100 cycles, and the statistics of the obtained V th distributions are reported in Fig. 8. We show almost no overlapping among SET and RESET distributions, without the use of smart programming or Program&Verify protocol.

Conclusions

A "Double-Patterned Self-Aligned" PCM-OTS structure was presented, targeting the co-integration of PCM "Wall" structure with OTS in the BEOL of the integration. The optimization of the fabrication process was discussed, presenting some of the implemented solutions. 1S1R devices, with dimensions down to 60 nm×80 nm, were realized and tested at statistical level in 1 kb 1T1S1R arrays. A reliable reading window of 1.5 V obtained with almost no SET/RESET distributions overlapping along 100 cycles.

Fig. 4 :

 4 Fig. 4: a) Resistance (median values) measured wrt Cross size on dedicated test structures (PCM layer w/o OTS i.e. 1R devices) for alignment control. b) Evolution of mean R V vs reciprocal of Wall size in 1R devices. The smallest Wall size of 40 nm shows a higher variability.

Fig. 2 :Fig. 5 :Fig. 7 :

 257 Fig.2: DPSA PCM-OTS process integration highlighting the "Wall" SA patterning (top) and the following "Cross" patterning (bottom).

SCU

Fig. 3 :Fig. 8 :Fig. 1 :Fig. 6 :

 3816 Fig.3: Contact recovery after "Wall" patterning and encapsulation: a) "Etch-back" approach: strong dependency on the devices density. In the TEM it is compared the nominal (good) density wrt 4x (bad). b) "Full CMP" approach: good uniformity, confirmed by SEM (section + BSE and SE).

Acknowledgments

This work was partially funded by European commission, French State and Auvergne-Rhone Alpes region through ECSEL-IA 101007321 project StorAIge and French Nano2022 program.


Philippe Blache
email: blache@ilcb.fr

Massina Abderrahmane
email: fm_abderrahmane@esi.dz

Stéphane Rauzy
email: stephane.rauzy@univ-amu.fr

Roxane Bertrand
email: roxane.bertrand@univ-amu.fr

An integrated model for predicting backchannel feedbacks

Keywords: Backchannels, embodied conversational agent, multimodal feedback, rule-based model

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Backchanneling constitutes one the keys for intelligent virtual agents naturalness, which is closely related to their capacity in generating prompt and adequate feedbacks to the speaker's production. Several predictive models have been proposed, involving different cues from prosody, gestures, lexicon, syntax or semantics. Unfortunately, acquiring all such cues in real time, which is mandatory in the case of conversational agents [START_REF] Kopp | An Architecture for Fluid Real-time Conversational Agents: Integrating Incremental Output Generation and Input Processing[END_REF], can be difficult. One simple solution consists in using temporal features [START_REF] Poppe | Backchannel strategies for artificial listeners[END_REF]. However, we also need to determine backchannels' type: generic or specific (initially called continuers or assessment) [START_REF] Bavelas | Listeners as co-narrators[END_REF] . These types involve different levels of processing, depending on whether they require or not deep understanding of the speaker's message. Some backchannels (esp. generic) are produced almost automatically and can be predicted by low-level cues. Some others require a certain level of semantic processing. The problem is that the production of these two types of backchannels is based on two different predictive mechanisms that are potentially in conflict because applied in parallel. We propose in this paper an approach avoiding this issue by implementing a single-route backchannel predictive model, generating appropriate backchannels in real time at a fine-grained level. The method we propose relies on a generic architecture. However, all dialogue systems being dependent from the task, we illustrate the approach with a specific application aiming at training human doctors to break bad news to virtual patients [START_REF] Ochs | Training doctors' social skills to break bad news: Evaluation of the impact of virtual environment displays on the sense of presence[END_REF]. In this system, the virtual patient is mainly a listener, capable of answering questions, requesting for clarification and producing backchannels.

BACKCHANNELS TYPES

Two BC categories are usually distinguished: generic (displaying attention to the speaker) and specific (expressing responses to the content of the speaker's production) [START_REF] Bavelas | Listeners as co-narrators[END_REF][START_REF] Bertrand | Co-narration in French conversation storytelling: A quantitative insight[END_REF][START_REF] Tolins | Addressee backchannels steer narrative development[END_REF]. They correspond to different communicative functions and can be expressed in different modalities: verbal, visual or multimodal. Our IVA being developed for French, we mainly focus the most frequent BCs given in [? ]: oui (yes) ah oui (oh yes), mh, d'accord (agree), ok, voilà (that's it), non (no), oh non (oh no), bon (well) ah bon (is that so?) (note that this list is comparable to that in English). In terms of types, oui and mh are typically generic whereas oh non or d'accord are specific.

Visual backchannels on their side correspond to many different types. They are mainly head movement (nod, jerk, shake, tilt, turn, waggle), facial expressions (smile, laughter) or eyebrow (frowning, raising). As noted in [START_REF] Ferré | Unimodal and Bimodal Backchannels in Conversational English[END_REF], they are less disruptive and in a large majority used as generic BCs. Types of backchannels Bimodal backchannels (involving both visual and verbal productions) are also very frequent and play an important role [START_REF] Ferré | Unimodal and Bimodal Backchannels in Conversational English[END_REF]. [START_REF] Bevacqua | Multimodal Backchannels for Embodied Conversational Agents[END_REF] has specifically studied a subset of such bimodals for IVA, associated with their functions: nod+yeah (agreement), shake+no (disagreement), smile+ok (interest), raise eyebrows+ooh (understanding), etc. In some cases, bimodality can reinforce the function: for example, bimodal BCs show a stronger agreement than unimodal ones. BC positions in the discourse also have an influence on their modality [START_REF] Ferré | Unimodal and Bimodal Backchannels in Conversational English[END_REF]: verbal BCs are preferably used during pauses, whereas visual BC are more likely to occur during speech. In terms of temporal features, verbal BCs are preferred soon after the beginning of a turn whereas visual BCs are usually produced later. We propose in this paper, beside the generic BC type, to precise the specific one by only using certain subtypes: agreement and disagreement (both being the most frequent in our corpus), completed with two other subtypes fitting with the specific purpose of our study: surprise and fear. The table 1 summarizes the different backchannels according to their modality and type.

BACKCHANNELS PREDICTIVE FEATURES

Many studies showed that the number of BCs has consequence on agent's naturalness. It is then important to identify as many BC triggers as possible for elaborating a natural model. In terms of temporal cues [START_REF] Poppe | Backchannel strategies for artificial listeners[END_REF] suggests that a BC may occur in average every 6.5 seconds. Another direct placement model consists in generating feedbacks systematically during breaks. These models are very simple, but produce unnatural behaviors of the agent. We have proposed more elaborated approach taking into account the temporal context in the framework of ACORFORMed doctor/patient dialogues [START_REF] Penteado | Evaluating Temporal Predictive Features for Virtual Patients Feedbacks[END_REF], based on the duration of the doctor's last silent pause, the duration since the doctor's last silent pause and the duration since last patient feedback. On their side, [START_REF] Poppe | Backchannel strategies for artificial listeners[END_REF][START_REF] Ward | Prosodic features which cue back-channel responses in English and Japanese[END_REF] have presented more detailed mechanisms based on prosodic and duration features:

• After a region of pitch less than the 26th-percentile pitch level, lasting 110ms, after 700ms of speech, with no BC within the preceding 800ms, after 700ms wait. • After a pause of 400ms preceded by at least 1000ms of speech, where the last 100ms contain a rising or falling pitch of at least 30Hz with no BC within the preceding 1400ms.

Several works [START_REF] Gravano | Turn-taking cues in task-oriented dialogue[END_REF][START_REF] Meena | Data-driven models for timing feedback responses in a Map Task dialogue system[END_REF] propose models based in particular on different prosodic backchannelling cues such as final rising intonation, higher intensity level, higher pitch level, inter-pausal duration, etc.

We propose to use a selection of such features (that can be recognized in real time): Prosody (time since IPU began, silent pause for over 200ms), Discourse (last dialogue act), Syntax (type of the last POS bigram), Semantic (semantic type of the last entity).

Moreover, we also identified sequences of multimodal features (mainly POS tags and gestures) that could be used as backchannel predicting cues [START_REF] Porhet | Mining a Multimodal Corpus of Doctor's Training for Virtual Patient's Feedbacks[END_REF]. From these sequences, several rules have been extracted for generating backchannels:

doctor_head_nod ⇒ patient_head_nod doctor_verb, doctor_head_nod, doctor_noun ⇒ patient_head_nod doctor_adverb ⇒ patient_head_nod doctor_medical_vocabulary ⇒ patient_head_nod 4 A PREDICTIVE RULE-BASED MODEL
Instead of a classical architecture based on two types of processing (an automatic loop generating low-level BC vs. a deeper one for semantic processing), we propose a unique model integrating both levels. This approach is not only more efficient, but also closer to a parallel human-like processing.

At the technical level, the different techniques now available in dialogue technology give access to many types of linguistic information in real time. More precisely, besides low-level classical cues (such as breaks, turn length, POS, etc.), we can also recognize during the dialogue semantic and discourse-level cues such as dialog acts, slot-based semantic structures as well as the position in the discourse ( [START_REF] Bertrand | Co-narration in French conversation storytelling: A quantitative insight[END_REF] shows how discourse phases influence BC types, generic BCs being favored when reaching the core of the narration). These different features all contribute to the generation of backchannels by deciding their type and placement. The interest of such an integrated model, on top of its homogeneity and its naturalness, is to avoid conflicts in the generation process, occurring when different backchannels can be generated simultaneously from two different loops.

Our proposal consists in integrating this BC generation model into an embodied conversational agent equipped with a dialogue system. The first question consists in identifying the span of each processing step that is traditionally a temporal window in the signal (e.g. each 3 seconds), an inter-pausal unit (i.e. between pauses of at least 200ms), a turn, etc. Embodied dialogue systems mainly have two input streams: the audio signal and its transcription (note that, in spite of the important role it can play, we put aside at this stage the gestural speaker's behavior in order to simplify the type of cues to be acquired in real time). The audio stream makes it possible to acquire temporal and prosodic features (silent pauses, pitch, IPU duration, etc.) where the transcription stream leads to the other linguistic cues. In this second stream, features from different domains can be acquired. First, at the lexical level, we have seen that morphosyntax (POS tags) plays an important role: for example, certain POS bigrams such as V-Adv favors the production of a generic BC. In the same vein, some terms (or more generally some semantic types) can be associated to a specific BC (surprise, fear, etc.). We also know that the introduction of new terms in the discourse (new referents) may be associate to a specific BC, indicating for example an agreement. At the discourse level, different studies have shown that the structure of the conversation and more precisely its phasing (opening, closing, argumentation, etc.) can also be associated with specific listener's reactions [START_REF] Bertrand | Co-narration in French conversation storytelling: A quantitative insight[END_REF]. This information can be directly obtained thanks to dialog act (DA) classification [START_REF] Chen | Dialogue Act Recognition via CRF-Attentive Structured Network[END_REF][START_REF] Stolcke | Dialogue act modeling for automatic tagging and recognition of conversational speech[END_REF] : in the case of our application domain, DAs correspond to dialogue phases. It is thus possible to trigger specific BCs upon phase change (e.g. when the doctor starts announcing a bad new). Finally, semantics plays a central role in generating specific BCs: many listener's reactions are triggered upon instantiation of the semantic structure, which is represented in task-oriented dialogues by a common ground (CG) [START_REF] Blache | Dialogue management in task-oriented dialogue systems[END_REF]. In such approaches, the CG consists in a set of structured frames made of different variables. The construction of the semantic structure consists in instantiating these variables during the discourse, which can elicit specific BCs (related to the content itself or to the global process of CG construction). We propose to implement a semantic-based BC generation by associating CG slots to specific BCs. For example, a BC expressing fear can be triggered when the slot "urgency" is instantiated.

The implementation of a unique loop for generating BC requires an appropriate segmentation of the input stream in order to take advantage of all different types of information we want to bring into the model. The processing span need to be large enough to capture semantic and discourse-level information, but not too long in order to allow real-time reactions. This means that neither arbitrary short segments nor entire IPUs returned by the speech recognition system (that can last in some cases more than 20 seconds) are adequate for our model. We propose an intermediate segmentation approach for answering these needs consisting in segmenting the input flow by using discourse markers (mais, donc, puis, alors, etc.) that indicate approximately a change between different discourse units. The list of such items being closed and very short, the mechanism simply consists in checking whether a marker appears. In such case, the current segment is parsed and the linguistic processing tools extract in real time the different features such as POS tags, semantic types of the substantives, dialogue acts and slot instantiation of the common ground in order to generate adequate BCs. Beside them, the different audio features are kept updated, in particular the duration since the last feedback, the indication of the current state of the production (speech or silent pause), the duration of the speech since the last pause, the duration of the pause, etc. Note that in the type of data we are working on, most of BCs are triggered by linguistic cues, only few of them being generated by temporal features. As a consequence, we give them a priority in the algorithm, temporal cues being only used as secondary triggers. Concretely, during the speech, we first look at linguistic BC cues whereas during a pause, the BC generator is mainly based on pause duration. In the case where none of these contexts trigger any BC, then the most general temporal cue (duration since the last BC) is applied. The algorithm 1 presents the general mechanism for controlling BC generation. This algorithm renders possible to process all types of BC generation in a single loop. It distinguishes different situations, hierarchically ordered (avoiding a concurrent processing generating possible conflicts). At the higher level, we have seen that slot filling (i.e. instantiation of the semantic structure) may directly triggers backchannels. In such case, the BC type associated with the slot becomes the current BC type value (e.g. agreement, surprise, etc.). If no BC comes from slot filling, then the different BC cues are extracted from the analysis of the current segment (e.g. dialogue acts, semantic type, POS sequences, etc.). These cues serve as input to a BC type identification function, based on a set of production rules. Note that the cues as well as the rules can be applied to the audio signal as well as the transcription (in other words, takes into account prosody, discourse, syntax etc.). Figure 2 gives examples of such rules. Given the BC type to be generated and the current mode (pause or speech), the last step consists in generating the BC itself, by choosing among a list of possible candidates (as described 1). Note that this list is in a probability space which also depends on the current state (visual or bimodal BC will be preferred during speech where verbal BCs will be favored during pauses).

CONCLUSION

We have presented in this paper a method for generating backchannels which is, at the difference with other approaches, integrated into a single loop. This approach takes advantage of different NLP technologies for a real time acquisition of different linguistic features involved into a high-level BC generating model. We have implemented this model into a virtual agent equipped with a dialogue system. The resulting application offers a very reactive artificial listener, capable of generating a large variety of multimodal BCs. This system is currently under evaluation.

Figure 1 :

 1 Figure1: Types of backchannels Bimodal backchannels (involving both visual and verbal productions) are also very frequent and play an important role[START_REF] Ferré | Unimodal and Bimodal Backchannels in Conversational English[END_REF].[START_REF] Bevacqua | Multimodal Backchannels for Embodied Conversational Agents[END_REF] has specifically studied a subset of such bimodals for IVA, associated with their functions: nod+yeah (agreement), shake+no (disagreement), smile+ok (interest), raise eyebrows+ooh (understanding), etc. In some cases, bimodality can reinforce the function: for example, bimodal BCs show a stronger agreement than unimodal ones.

Figure 2 :

 2 Figure 2: BC generation rules


Samir Beneddine

Nonlinear input feature reduction for data-based physical modeling

Keywords: deep learning, mutual information, dimensional analysis, physical modeling

This work introduces a novel methodology to derive physical scalings for input features from data. The approach developed in this article relies on the maximization of mutual information to derive optimal nonlinear combinations of input features. These combinations are both adapted to physicsrelated models and interpretable (in a symbolic way). The algorithm is presented in detail, then tested on a synthetic toy model. The results show that our approach can effectively construct relevant combinations by analyzing a strongly noisy nonlinear dataset. These results are promising and may significantly help training data-driven models. Finally, the last part of the paper introduces a way colorblackto account for the physical dimension of data. The test case is a synthetic dataset inspired by the Law of the Wall from turbulent boundary layer theory. Once again, the algorithm shows that it can recover relevant nondimensional variables colorblackfor data-base modeling.

Introduction

Open physical modeling questions have recently benefited from the accelerating developments of Machine Learning tools. In particular, Deep Learning (DL) carries high hopes of tackling numerous unresolved physical problems. For instance, for the specific field of fluid mechanics, the work of [START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF] has been among the precursors for many papers attempting to propose new DL-augmented Reynolds-Averaged Navier-Stokes (RANS) models [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework[END_REF][START_REF] Volpiani | Machine learning-augmented turbulence modeling for rans simulations of massively separated flows[END_REF]. To mention a few other applications, DL has also been considered as a tool to produce new Wall Models [START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF], Subgrid-Scale Models [START_REF] Vollant | Subgrid-scale scalar flux modelling based on optimal estimation theory and machine-learning procedures[END_REF], and Transition Models [START_REF] Yang | Improving the k-ω-γ-ar transition model by the field inversion and machine learning framework[END_REF]. Needless to say that this scientific trend goes well beyond the sole topic of fluid mechanics or even physics and concerns virtually all open modeling problems (for instance, in biotechnology [START_REF] Gao | Deep learning in protein structural modeling and design[END_REF], solid mechanics [START_REF] Haghighat | A physicsinformed deep learning framework for inversion and surrogate modeling in solid mechanics[END_REF], molecular dynamics [START_REF] Zhang | Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics[END_REF], epidemiology [START_REF] Shorten | Deep learning applications for covid-19[END_REF], etc.).

Yet, relying on a purely data-based approach has shown mitigated results. To our knowledge, no data-driven model has definitively answered any of the topics mentioned above. Indeed, neural networks (NN) are known to have poor extrapolation capabilities [START_REF] Hettiarachchi | The extrapolation of artificial neural networks for the modelling of rainfall-runoff relationships[END_REF], which results in low accuracy when a DL-model is used for physical conditions that have been unseen during training. Additionally, NNs are often black boxes from which it is hard to gain new scientific knowledge. This lack of generality and understandability of DL models may be one of the core challenges in modern Machine Learning applied to physical problems.

Consequently, recent papers have attempted to incorporate physical knowledge into these data-based approaches. Several works have directly included the governing equations of a dynamical system within the loss function of neural networks. This led for instance to the so-called Physics-Informed Neural Networks (PINN) [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF], used for a wide range of physics-related problems [START_REF] Cai | Physics-informed neural networks (pinns) for fluid mechanics: A review[END_REF][START_REF] Mao | Physics-informed neural networks for high-speed flows[END_REF][START_REF] Sahli Costabal | Physics-informed neural networks for cardiac activation mapping[END_REF]. Another way to include physics within machine learning approaches has been done through the input features. For instance, [START_REF] Yang | Predictive large-eddy-simulation wall modeling via physics-informed neural networks[END_REF] showed in their works that scaling input features based on prior knowledge of the physics involved leads to far better results than just blindly feeding the NN with all available quantities. Similarly, [START_REF] Volpiani | Machine learning-augmented turbulence modeling for rans simulations of massively separated flows[END_REF] wrote in their paper that "the choice of the input features is crucial" to obtain data-based RANS models with acceptable fidelity. These are just two examples of many works demonstrating that successful DL approaches for physical modeling depend on the judicious choice of input quantities. Therefore, the adequate definition of the input features is a central question for data-driven modeling.

Note that feature selection is a research topic highly active in the Machine Learning community, especially on non-physics-related topics and classification tasks, with plenty of articles dedicated to the question (to name a few, [START_REF] Dash | Feature selection for classification[END_REF][START_REF] Blum | Selection of relevant features and examples in machine learning[END_REF][START_REF] Battiti | Using mutual information for selecting features in supervised neural net learning[END_REF][START_REF] Sindhwani | Feature selection in mlps and svms based on maximum output information[END_REF][START_REF] Bollacker | Linear feature extractors based on mutual information[END_REF][START_REF] Tadist | Feature selection methods and genomic big data: a systematic review[END_REF]). Researchers have gathered existing methods into five main categories: filter methods, wrapper methods, embedded methods, ensemble methods, and integrative methods (see, for instance, [START_REF] Naik | A novel sensitivity-based method for feature selection[END_REF][START_REF] Tadist | Feature selection methods and genomic big data: a systematic review[END_REF][START_REF] Bolón-Canedo | A review of feature selection methods on synthetic data[END_REF] for details). For this work, we are particularly interested in finding methods to identify a priori (before any training) some relevant features just by analyzing existing data. The techniques developed in this article can be viewed as assistive tools for classical or data-driven modeling. It falls into the filter methods category, where numerous approaches have been considered based on linear correlation, Fisher score, etc. (see, for instance, [START_REF] Naik | A novel sensitivity-based method for feature selection[END_REF]).

Several particularly promising approaches for filter methods rely on information theory. Mutual information maximization (a notion explained in the paper) has been extensively used to eliminate redundant inputs among large sets of inputs. For instance, information gain is a standard tool for text classification [START_REF] Yang | A comparative study on feature selection in text categorization[END_REF], where mutual information between a single feature and a class variable estimates the relevance of each feature one by one. Along the same line, [START_REF] Battiti | Using mutual information for selecting features in supervised neural net learning[END_REF] developed the mutual information feature selection (MIFS) algorithm for classification problems that searches greedily a set of features with high-mutual information with class labels but low mutual information among chosen features. colorblack Other similar algorithms exist, such as the Maximum Relevance Minimum Redundancy (MRMR) approach proposed by [START_REF] Peng | Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy[END_REF], or the Joint Mutual Information (JMI) method from [START_REF] Yang | Feature selection based on joint mutual information[END_REF] to derive relevant subsets of features. These approaches allow for the elimination of redundant information within features. Among other works, one may cite [START_REF] Bollacker | Linear feature extractors based on mutual information[END_REF], which used mutual information between class variables to form optimal linear combinations of inputs. However, as mentioned by [START_REF] Sindhwani | Feature selection in mlps and svms based on maximum output information[END_REF], these methods, which deal with continuous variables in classification tasks, demand large amounts of data and high-computational complexity. colorblackMany other algorithms exist and they have been widely used in particular for classification tasks (e.g. medical image analysis [START_REF] Allam | A study on optimization techniques in feature selection for medical image analysis[END_REF], phenotype classification [START_REF] Ding | Minimum redundancy feature selection from microarray gene expression data[END_REF], etc.). A complete review may be found in [START_REF] Li | Feature selection: A data perspective[END_REF].

The present work is inspired by this existing literature, except that the developed approach is tailored explicitly for physics-related problems. In contrast to the articles cited above, this study focuses on regression tasks, and we do not want to perform a simple elimination of input features. Instead, it aims to form specific nonlinear feature combinations that are both relevant for physical modeling and understandable (in a symbolic way). Additionally, the last part of the paper proposes an approach that accounts for the physical dimension of input features. The introduced algorithms are new, and they come with several implementation challenges that are addressed in the paper. Finally, one major novelty is that all methods are based on the recent work of [START_REF] Belghazi | Mutual information neural estimation[END_REF], which provides innovative techniques to estimate the mutual information between two random variables.

The paper is organized as follows. The first part defines the context of the study and illustrates the issues addressed in this article using a synthetic toy model. The complete approach and corresponding algorithm are then detailed. This algorithm is tested in the next section on the previously introduced toy model. A variant of this model is also considered to present a methodology that produces several reduced variables instead of a single one. Finally, before concluding, the last section presents a way to colorblack include dimensional knowledge about the input features in the algorithm. The method is tested on a synthetic case inspired by the Law of the Wall from the turbulent boundary layer theory. All results of the paper have been produced using the pyTorch library [START_REF] Ketkar | Introduction to pytorch[END_REF].

Definitions and proposed strategy

Notations and definitions

Let us consider a quantity f = (f 0 , f 1 , . . . , f n ) ∈ R n to model. The set of input variables of the sought model is q = (q 0 , q 1 , . . . , q m ) ∈ R m . In the context of fluid dynamics, f may be for instance a corrective term for a RANS model (as done in [START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF] or [START_REF] Volpiani | Machine learning-augmented turbulence modeling for rans simulations of massively separated flows[END_REF]), and q may be some local fluid variables such as the density, velocity, molecular viscosity, shear stress tensor components, etc.

Let us assume that f and q are linked through an unknown stochastic model M: M(q, η) = f , with η a stochastic variable associated with some incompressible noise (that is not modeled). The initial variable set is assumed to be not optimal, i.e., it is hard to learn a data-based model M from (q 0 , q 1 , . . . , q m ). It may be because this set of inputs contains redundant or useless information or because there exists a reduced set of variables that makes the mathematical model simpler (this aspect is further developed in section 2.2). Such improper input quantities will be called "naive" input variables in the following. As mentioned in the introduction, choosing physically relevant input features is a common issue for data-driven physical modeling. More generally, the proper selection of inputs is a known problem in deep learning that may drastically affect the convergence, generality, and accuracy of the learned model (see section 2.2 as an illustration of this).

Finally, let us say that some reduced sets of more relevant input features exist. In physics, this would typically be a set of non-dimensional quantities or dimensional quantities involving an appropriate scaling. Therefore, the unknown model M may be split into two functions E and M ′ , such that M = M ′ •E, with E : R n → R l the mapping of the naive input vector q to a bettersuited l-dimensional space that feed the model M ′ . This new model M ′ is ideally easier to learn from data and may have enhanced understandability. The black square represents the range of values for (q 0 , q 1 ) used for training neural networks to estimate f (q 1 = 0 excluded).

Illustrative example 2.2.1. Definition of the toy model

This section presents a toy model used for a regression task using deep learning. It illustrates some intuitive ideas mentioned earlier, and it will be used later in the article to test the new approaches developed in this work.

Let us consider a database containing 4000 triplets (q

(i) 0 , q (i)
1 , f (i) ) obtained from the following nonlinear stochastic toy model f = M(q 0 , q 1 , η) = q 2 0 q 1 + 3 cos(2π

q 2 0 q 1 )(1 + η), (1)
with η a white noise of amplitude 0.5, and (q

(i) 0 , q (i) 1 ) taking values in [-1, 1]× [-1, 1] \ {0}.
A graphic representation of this function is shown in figure 1. This toy model has been designed to be strongly nonlinear, very noisy, and such that a simple regression task by a neural network may become challenging if not done adequately. This model is particularly interesting due to its behavior near q 1 = 0: f diverges and oscillates very rapidly in this neighborhood, causing difficulties to perform a regression using a neural network.

A first "naive" strategy to predict f consists of training a neural network N directly from (q 0 , q 1 ) using available data. Alternatively, since q = q 2 0 q 1 is a reduced variable that changes M into a simpler model M ′ :

f = M ′ (q, η) = (q + 3) cos(2π q)(1 + η), (2)
one may train a network Ñ to estimate f from q using the same training data. This approach is called hereafter a "model-informed" strategy (since it requires some prior knowledge of the model). We will demonstrate some clear advantages of this latter strategy in the following.

Comparison of the naive and model-informed regression

To perform a fair comparison between the "naive" and the "model-informed" strategy, N and Ñ have the same number of hidden layers, units per layer, and activation functions. Training is performed the same way: the database is split in two for the validation (20% of samples) and training (80% of samples), and an Adam optimizer [START_REF] Kingma | A method for stochastic optimization[END_REF] is used with the same learning rate to minimize the mean square error between the network output and the actual value of f . No extra penalization or more advanced technique is used here. Training is performed until an early-stopping criteria (based on the validation loss) is met. Full details are given in Appendix A.

Learning curves from figure 2 show that the naive strategy has a poor and slow convergence. By comparison, the model-informed approach is much easier to train. The interpolation and extrapolation capability of each network may be seen in figure 3. The naive approach is not only unable to extrapolate outside from the training range (q

(i) 0 , q (i) 1 ) ∈ [-1, 1] × [-1, 1] \ {0}
, but even its interpolation capabilities are unsatisfactory. On the other hand, the model-informed strategy has much better interpolation capabilities (figure 3(right)). It is even able to extrapolate to some extent: while the center-left and center-right parts of the domain in figure 3(right) are not well predicted, the rest shows good accuracy even for unseen values of (q 0 , q 1 ), because the corresponding values of q have actually been seen during training.

These results do not demonstrate that one cannot develop a more advanced training strategy or network architecture to make the naive strategy work with this database. Still, they highlight that using specific input features may significantly ease training. This numerical experiment is a simple illustration of known neural network training behavior. Yet, it highlights well the motivations behind the present paper: the need for a strategy to deduce the "proper" reduced variable(s) to use from a database. For instance, for this toy model, it may ideally tell to use q 2 0 q 1 as an input feature. This would ease the network's training and improve its interpolation/extrapolation capabilities. -2 The black square represents the range of values for (q 0 , q 1 ) used for training (everything outside the square shows extrapolation capabilities).

0 2 q 0 -2 -1 0 1 2 q 1 -5 0 5 -2 0 2 q 0 -2 -1 0 1 2 q 1 -5 0 5

7

-2

0 2 q 0 -2 -1 0 1 2 q 1 -5 0 5 -2 0 2 q 0 -2 -1 0 1 2 q 1 -5 0 5
Figure 4: (Left): estimation of f obtained using by q-1 as input feature. (Right): estimation of f obtained using by q2 as input feature. The black square represents the range of values for (q 0 , q 1 ) used for training (everything outside the square shows extrapolation capabilities).

Non-unicity of reduced variables

The previous results show that the reduced variable q is a much better input feature choice than (q 0 , q 1 ). But one may wonder if other choices could be even better. In particular, for this paper, it is interesting to see if variables of the form qx are also good choices (the reason for this particular form is explained in section 3.2).

In general, as soon as x ̸ = 0, this would make a better input feature than (q 0 , q 1 ): the dimensional reduction of the input space mechanically gives a network trained with such a variable extrapolation capabilities to some extent. The only issue may be that for non-integer values of x, qx is not be defined for negative values of q 0 or q 1 , so one may favor integer values for x.

Nonetheless, most variables of the form qx yield nearly equivalent results in terms of interpolation / extrapolation performances. For instance, figure 4(left) shows the result field obtained from a network trained with the input feature q-1 : except for the zone near q 0 = 0 (harder to learn due to the discontinuous behavior induced by this input feature choice), the overall result is similar to that from the previous section. In contrast, using q2 provides significantly better results (figure 4(right)). This is because the relation between (q 0 , q 1 ) and f displays some symmetries, but the network is unaware of this and has to learn it from data. When using q2 as an input feature, the reduced model to learn M ′ has no longer any symmetry. Thus the network has less information to learn and reaches better regression capabilities.

The conclusion will further discuss these few remarks since they pro-vide interesting insights and maybe future work directions about the general question of choosing proper input features, which is the central topic of this article.

Remark on correlation analysis to find relevant input features

One widely used tool to determine the relevancy of a given input variable to model another is the computation of correlation coefficients. The choice of this particular toy model is interesting because it shows that a correlation analysis cannot indicate that q is a well-suited variable: even if η = 0, the Pearson correlation coefficient between q and f for (q

(i) 0 , q (i) 1 ) ∈ [-1, 1] × [-1, 1] \ {0}
would only be 0.04. This limitation is due to the linear nature of a correlation analysis, which would lead here to discard q as a relevant input variable to model f .

Mutual information concept

Before introducing a new approach for finding good input features for physical models, let us define the notion of "proper" input variable better. This paper bases the answer on information theory: a well-chosen input feature should provide as much information as possible on the quantity to model. This can be quantified using mutual information, a concept defined below.

In information theory, the mutual information I(X, Y ) between two random variables X, Y from a space X × Y is defined as

I(X, Y ) = D KL (P Y,X ∥P X ⊗ P Y ) , (3)
where P (X,Y ) is the joint probability distribution of the pair (X, Y ), P X and P Y are the marginal distribution of the pair (X, Y ), and D KL (•∥•) is the Kullback-Leibler (KL) divergence. After expending this expression into differential form

I(X, Y ) = X ,Y P X (X)P Y |X (Y ) log P Y |X (Y ) P Y (Y ) dXdY = - Y P Y (X) log P Y (Y )dY + X P X (X) Y P Y |X (Y ) log P Y |X (Y )dY dX = E Y ∼P Y [-log P Y (Y )] -E X∼P X ,Y ∼P Y |X [-log P Y |X (Y )], (4)
where P Y |X is the conditional distribution of Y knowing X, a well-known alternative definition of I involving Shannon's entropy H appears

I(X, Y ) = H(Y ) -H(Y |X). (5)
Equation ( 5) gives an alternative view on I: it is linked to the remaining uncertainty on Y (respectively X) once X (respectively Y ) is known. In other words, I(X, Y ) is the amount of information contained in one random variable about the other. For instance, if the two variables are independent, then H(Y |X) = H(Y ) and I(X, Y ) = 0. Contrary to linear correlation, it is a measure of true dependence since it is able to correctly quantify nonlinear statistical dependency [START_REF] Kinney | Equitability, mutual information, and the maximal information coefficient[END_REF] (which is not the case for linear correlation as illustrated with the toy problem from section 2.2).

The proposed strategy of this paper, detailed in the next section, relies on maximizing the mutual information between some input combinations and the quantity to model: some parameters A (the weights of a neural network E A ) will be optimized to maximize I(E A (q), f ) (using notations from 2.1), thus providing the mapping from naive to relevant input features.

Mutual information maximization strategy

In most cases, it is impossible to directly maximize I(E A (q), f ) because it involves the unknown posterior distribution P q|f (using the notation q = E A (q)). To overcome this issue, we propose an approach based on the recent work of [START_REF] Belghazi | Mutual information neural estimation[END_REF], where the mutual information is first estimated using a dual representation of the KL-Divergence (providing a lower bound). This estimation is then used to maximize the mutual information.

The estimation of the mutual information relies on the Donsker-Varadhan (DV) representation [START_REF] Donsker | Asymptotic evaluation of certain markov process expectations for large time[END_REF], based on the following theorem (see for instance [START_REF] Belghazi | Mutual information neural estimation[END_REF] for the proof)

Theorem 1.
Let Ω be a sample space, and P and Q two given probability distributions on Ω. Then, the KL-divergence admits the following dual representation:

D KL (P ∥Q) = sup T :Ω→R E P [T ] -log(E Q [e T ]), (6)
with the supremum taken over all functions such that the expected values are finite.

Since equation (3) defines the mutual information as the KL-divergence of the joint distribution and the product of the marginals, the DV-representation straightforwardly provides lower bounds that may be used to maximize the mutual information I(E A (q), f ).

To estimate I(X, Y ), the idea is to consider a large family of function T ψ : X , Y → R parametrized as a neural network with parameters ψ ∈ R K (with K the number of weights and bias of the network), and set ψ such that it maximizes the quantity

Ĩψ (X, Y ) = E P X,Y [T ψ ] -log(E P X ⊗P Y [e T ψ ]). ( 7
)
This quantity is a lower bound for I(X, Y ), and such a network is called a statistics network in the Mutual Information Neural Estimator (MINE) framework developed by [START_REF] Belghazi | Mutual information neural estimation[END_REF]. Given the universal approximation theorem for neural networks [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], one may expect to get through Ĩψ (X, Y ) an arbitrarily tight bound given that the network complexity is high enough. Details about the network architecture for T ψ used in the paper are given in the appendices.

Adequate choice of function space to represent input features

This section proposes a particular network architecture for E A suited for physical modeling, and that provides interpretable results. When looking at existing physical models, relevant input variables come from scalings of the form q = i q α i i .

It generally does not involve more complex functions for dimensional reasons.

For instance, the logarithmic Law of the Wall, which models a part of the velocity profile of turbulent boundary layers in fluid mechanics, is a relation linking the streamwise velocity u of the flow to the distance from the wall y that reads

u + = 1 κ ln(y + ) + C + , (9)
with

τ w = µ ∂u ∂y y=0 , u * = τ w ρ , y + = yρu * µ , u + = u u *
where ρ is the fluid density, µ the dynamic viscosity, and κ and C + two constant values. This law provides u from the intermediary variables y + and u * , which are combinations of the simpler variables y, ρ, µ, ∂u ∂y y=0 of the form [START_REF] Haghighat | A physicsinformed deep learning framework for inversion and surrogate modeling in solid mechanics[END_REF].

Therefore, we propose to use a logarithm representation of the quantity to model (we focus on log(f ) instead of f ), such that the input features are searched under the form

q = i α i log(q i ), (10)
instead of using equation [START_REF] Haghighat | A physicsinformed deep learning framework for inversion and surrogate modeling in solid mechanics[END_REF]. A neural network forming such quantities can easily be designed: it consists of a log activation followed by a linear layer with no bias. The weights of such a network are directly the exponents of equation ( 8). This architecture, called in the following a logarithmic network (LN), is represented in figure 5. Note that the inputs are first made positive by taking their absolute value to handle negative numbers properly. This has no impact in the method since it is the same set of exponents a i that relates q with q i and |q| with |q i |. As a side remark, an alternative solution could have been to use the complex definition of the logarithmic function

log(z = |z|e iθ ) = log(|z|) + iθ. (11)
One limitation remains: the LN cannot handle null values, which have to be removed or replaced by some ϵ value in the training database considered. Note that although it is not explored in this paper, by adding a linear layer before the log activation, this approach may be straightforwardly extended to the more general form of input features which may encompass some specific cases of physical scaling that may be not covered by relation [START_REF] Haghighat | A physicsinformed deep learning framework for inversion and surrogate modeling in solid mechanics[END_REF]. In any case, different functional spaces that would be more fitted to a given problem may easily be designed and used in the framework presented in this paper. Note that in the following, LN with a single output will be exclusively used (corresponding to m = 1 in figure 5) due to the particular strategy used for models with multiple reduced variables (explained in section 3.3).

q = i j a ij q j α i , (12)
q 0 q 1 q n . . . log(|q 0 |) log(|q 1 |) log(|q n |) . . . log(| • |) Linear layer n × m matrix A n i=1 a 0i log(|q i |) n i=1 a 1i log(|q i |)

Feature Discovery algorithm

This section combines the notion and techniques introduced above to propose the Feature Discovery (FeDis) algorithm. The idea is to maximize the mutual information between the output of an LN and a quantity to model using a MINE statistical network. The overall procedure to compute a single reduced variable is presented in algorithm 1. One may see that a penalization of the LN weights is added to the loss function. In the cases treated in the paper, a L 1 penalization was systematically considered since it promotes sparsity (and therefore helps to discard naive inputs that do not contain information about the modeled quantity) and was found to improve convergence overall. Other more advanced regularizations to handle multiple reduced variables are discussed in section 3.3.

Additionally, algorithm 1 shows one more implementation specificity: instead of maximizing the mutual information between E(q) and f , the algorithm uses E(q) and log(|f |). This has no impact on the algorithm (maximizing the mutual information between these transformed variables still provides the input that contains the most information about f ) and has the advantage Algorithm 1 Feature Discovery (FeDis) algorithm m ← 1 Set the dimension of the latent space for the reduced variables to 1 A, ψ ← Initialize network parameters for the 2 networks E A (LN network, output dim. = m) and T ψ while loss L not converged do Draw a batch of randomly sampled pair (q 1 , f 1 ), . . . , (q b , f b ) Eliminates/process pairs with (q i ) having null component(s) Form the pairs

c 1 = (E A (q 1 ), log(|f 1 |)), . . . , c b = (E A (q b ), log(|f b |)) (joint distribution) Form (f ′ 1 , . . . , f ′ b ) by shuffling (f 1 , . . . , f b ) (Marginal distribution) Form the pairs c ′ 1 = (E A (q 1 ), log(|f ′ 1 |))), . . . , c ′ b = (E A (q b ), log(|f ′ b |)) Evaluate the cost function L ← -1 b b i=1 T ψ (c i ) + log( 1 b b i=1 e T ψ (c ′ i ) ) [Optional] Add regularization terms L ← L + R(θ)
Jointly update θ and ψ to minimize L (gradient-based optimization) end while of yielding entries for T ψ of lesser magnitude (which helps avoid float-overflow issues that may occur when evaluating e T ψ (•) ).

Note that the original paper of [START_REF] Belghazi | Mutual information neural estimation[END_REF] about MINE networks introduced some advanced training techniques such as gradient clipping and a modified loss formulation to avoid bias during the gradient descent. These techniques have not brought any significant improvements to the results of this paper. They are therefore neither introduced nor used for the present article.

Test of the FeDis algorithm

Results on the augmented toy model

The toy model from section 2.2 is used to demonstrate the ability of the FeDis approach to finding relevant reduced inputs from data generated by a noisy nonlinear model. We consider an augmented set of 14 naive input features (q 0 , q 1 , . . . , q 13 ) that add 12 useless inputs to the original dataset (to assess the ability of the approach to sort useful/useless variables). The model is therefore f = M(q 0 , q 1 , q 2 , . . . , q 13 , η) = q 2 0 q 1 + 3 cos(2π

q 2 0 q 1 )(1 + η), (13)
with η a white noise of amplitude 0.5. Consistently with section 2.2, the training database is made of 4000 samples (q (i) 0 , . . . , q (i) [START_REF] Cai | Physics-informed neural networks (pinns) for fluid mechanics: A review[END_REF] , f (i) ), with each input q j randomly sampled in [-1, 1] \ {0}. All details to reproduce the test case are given in the appendices.

Figure 6 shows the corresponding evolution of the exponent value of each input during training. As expected, all useless variables are quickly discarded, and the algorithm yields the reduced variable q 3.01 1 /q 5.99 0 (i.e., approximately q-3 ). 0 200 400 600 800 epoch q 0 q 1 q 2 q 3 q 4 q 5 q 6 q 7 q 8 q 9 q 10 q 11 q 12 q 13 -5.99 +3.01 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 final values

Normalization of the exponents

The FeDis approach does not yield q = q 2 0 /q 1 , but instead a variable of the form qx . This was expected since I(q, f ) = I(q x , f ) for x ̸ = 0 (the mutual information I(X, Y ) is preserved by any invertible deterministic transformation of X, see [START_REF] Belghazi | Mutual information neural estimation[END_REF]). Therefore, given the structure of the LN network, the maximization strategy followed here leaves the exponent x undetermined (the convergence of x is eventually due to the L 1 penalization of weights that forbid the exponents to become too large). Finding exactly q 2 0 /q 1 would have been coincidental. But as mentioned in section 2.2, using qx as input feature provides the same advantages as using q. The only limitation is that non-integer exponents are not defined for negative numbers. But having this in mind, it is easy to rescale a posteriori the exponent to overcome this issue since the method provides a symbolic representation of the input. One could also imagine penalizing non-integer exponents (for instance, using the function p(a) = -λ int cos(2πa)), but that would add one extra hyperparameter λ int for an uncertain and arguable improvement of the original algorithm.

A simple a posteriori normalization to remove the undetermined nature of the exponent is given by algorithm 2, which sets the smallest exponents to 1 (discarding nearly zero values). If applied to the results obtained in section 3.1, the resulting reduced variable is q-1 = q 1 /q 2 0 .

Algorithm 2 Exponent normalization procedure NormalizeExponents(a = (a 0 , a 1 , . . . , a n-1 ))

a max ← max i (|a i |)

Dealing with multiple reduced variables

Only one reduced variable was needed for the toy model of section 3.1. But more than a single reduced variable may be needed for an actual physical model. To address this question, let us consider the following modified toy model f = M(q 0 , q 1 , q 2 , . . . , q 13 , η) = q 2 0 q 1 + 3 cos(2πq 0 q 1 )(1 + η).

Then, one may want to use reduced variables of the form (q a 0 , qb 1 ), with q0 = q 0 q 1 and q1 = q 2 0 q 1 . A first idea would be to use an LN network with an output dimension m = 2 (see section 2.5) such that the FeDis algorithm would yield a 2-dimensional vector whose components are feature combinations maximizing the information about f . In practice, this idea does not work. Indeed, the most straightforward pair of variables maximizing I(•, f ) is the naive couple (q 0 , q 1 ) (the deterministic part of f is fully defined once these two variables are known). Therefore, since the algorithm promotes sparsity, the result is likely to be (q 0 , q 1 ) (or a trivial variant). Actually, the algorithm could output nearly any pair of independent combinations of q 0 and q 1 and would still maximize the mutual information with f . We need a way of promoting a combination set where each component taken individually also maximize the mutual information.

The solution is to compute one reduced variable at a time. A first run of the algorithm provides the nonlinear combination of features that maximizes I(•, f ) (standard FeDis procedure). Then, the algorithm is rerun to give a different nonlinear combination that maximizes the mutual information, and so on. This one-by-one procedure forbids the production of trivial feature combinations such as (q 0 , q 1 ) and ensures that at each step, the next reduced variable is the best nonlinear combination given the previous ones. For this approach to work, at each step, an additional penalization term needs to be added to forbid the algorithm from producing over and over the same combination of features.

To illustrate the proposed approach, let us consider the second toy model defined by equation [START_REF] Mao | Physics-informed neural networks for high-speed flows[END_REF]. We proceed the same way as before: the training database is made of 4000 samples (q (i) 0 , . . . , q (i) [START_REF] Cai | Physics-informed neural networks (pinns) for fluid mechanics: A review[END_REF] , f (i) ), with each input q j randomly sampled in [-1, 1] \ {0}. The database is completed by the corresponding values of f generated from equation ( 14) with η = 0.5. Then, reusing the same techniques as before, the FeDis algorithm is run to get a first reduced variable (details regarding the neural networks are given in Appendix C). The result is shown in figure 7(top). The algorithm yields the reduced variable q -3.9 0 q -4.1 1 , which after normalization (defined in section 3.2) gives q0 ≈ q 0 q 1 (note that if rerun and initialized differently, the network would sometimes converged toward the second reduced variable (

q 2 0 q 1 ) x ).
Then, the algorithm is rerun with a regularization promoting different feature combinations. Note that any set of exponents that would be a multiple of those already obtained are unwanted (the exponents maximizing the mutual information are defined up to a multiplicative constant). The additional penalization term to achieve this is the following. Consider that the set of exponents obtained during the first run of the algorithm is the vector a = (a 0 , a 1 , . . . , a n ). Proportionality between the current exponents A = (A 0 , A 1 , . . . , A n ) during the rerun and a may be quantified using the normalized scalar product s s = ⟨a, A⟩ ∥a∥∥A∥

with ⟨•, •⟩ the classical dot product, and ∥ • ∥ the associated norm. When a and A are orthogonal, s = 0, and when they are aligned, s = 1. Then, the following penalization is added to the loss function of the FeDis algorithm

L a = λ a e -(s-1) 2 σ 2 , ( 16
)
0 200 400 600 800 1000 1200 1400 epoch q q q q q q q q q q q q q q -3.89 -4.10 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 +0.00 final values epoch q q q q q q q q q q q q q q -2.78 +1.40 +0.00 -0.00 +0.00 -0.00 -0.00 +0.00 -0.00 +0.00 -0.00 -0.00 -0.00 -0.00 final values with λ a a weighting factor for the penalization. Note that the use of a Gaussian function instead of s directly allows to not promote orthogonality of a and A (which is unwanted), but to only penalize too significant alignments. For the present study, σ 2 = 0.2. An exhaustive study on the optimal choice for σ 2 has not be conducted, but other values near 0.2 (0.15, 0.25) have been tested and they nearly did not changed the results.

The results obtained using this penalization are shown in figure 7(bottom) (implementation details given in Appendix C). The algorithm behaves as expected, and yields the second reduced variable q -2.8 0 q 1. 4 1 , which gives after normalization q0 = q 1 /q 2 0 .

Dimension-aware algorithm

Finding relevant reduced variables to produce a physical model echoes the well-studied question of dimensional analysis. The Buckingham π-theorem states that any physically meaningful governing equation involving k input variables can be reduced to an equation involving kl dimensionless parameters, with l the number of physical dimensions involved. Therefore, it is natural to see if the techniques introduced in this paper can be used to automatically perform a feature reduction accounting for the physical dimension of the input variables.

Physics-inspired test case

The following proposes an algorithm that answers the question mentioned above. To introduce it, let us focus on the already-presented Law of the Wall, a relation between the streamwise velocity u in turbulent boundary layers and the distance from the wall y, which links two dimensionless quantities y + and u + defined as

y + = y √ ρ ∂u ∂y y=0 √ µ , (17)
u + = u √ ρ µ ∂u ∂y y=0 , ( 18
)
where ρ is the fluid density, µ the dynamic viscosity, and κ and C + two constant values. The relation, already detailed in section 2.5, reads

u + = 1 κ ln(y + ) + C + . (19)
ω 0 µ LN network Update A 1 , . . . , A 7 (gradient-based minimization) Statistical network T Ψ X y = ρ ω 3 ω 2 ω 1 u ∂u ∂y (0) ω 0 µ X u = ρ ω 3 ω 2 ω 1 1 A 1 A 4 A 2 A 3 A 7 A 6 A 5 f y = log(|y|) + 7 i=1 A i log(|X y,i |) LN network 1 B 1 B 4 B 2 B 3 B 7 B 6 B 5 f u = log(|u|) + 7 i=1 B i log(|X u,i |) Loss Evaluation L = -ĨΨ (f y , f u )+ Regularization Update B 1 , . . . , B 7 (gradient-based minimization)
Update Ψ (gradient-based minimization)

Figure 8: Schematic representation of the FeDis variant algorithm that accounts for physical dimensions, applied to the synthetic Law of the Wall case.

In the following, we consider synthetic data generated using this relation. A database made of 10000 set of values (y, ∂u ∂y y=0 , µ, ρ, ω 0 , ω 1 , ω 2 , ω 3 ) randomly sampled from ]0, 1] 8 has been generated. Note that this range of values has been chosen for simplicity since the case serves illustrative and demonstration purposes (thus the name of "physics-inspired" test case). It does not correspond to realistic values for each feature. Then, u is evaluated for each sample using the following relation

u = u * 0.41 ln(y + ) + 5 (1 + η), (20)
with η a random white noise of amplitude 0.2. The chosen values κ = 0.41 and C + = 5 are the standard constant values. Note that the variables ω i are extra dummy quantities that are unused to compute u (they allow to assess the ability of the approach to discard useless data).

Length L Time T Mass M y 1 0 0

∂u ∂y 0 0 -1 0 µ -1 -1 1 ρ -3 0 1 ω 0 1 0 0 ω 1 0 0 1 ω 2 0 -1 0 ω 3 1 1 1 Length L Time T Mass M u 1 -1 0 ∂u ∂y 0 0 -1 0 µ -1 -1 1 ρ -3 0 1 ω 0 1 0 0 ω 1 0 0 1 ω 2 0 -1 0 ω 3 1 1 1
Table 1: Dimension of the features defining q y and q u . These tables define the so-called dimension matrices M d,y (left) and M d,u (right) associated with q y and q u respectively. They are used to promote dimensionless combinations.

is unwanted because it overlooks the relation between u and y. Therefore, the alignment of à and B needs to be penalized. This is done using the same alignment penalization from section 3.3 (equation ( 16)).

Results

Details and hyperparameters used to produce the results are given in Appendix D. Figure 9 shows the evolution of the exponents A and B during training. One may see that it produces the wanted result: the algorithm is able to nearly recover the expressions of y + and u + (equations ( 17) and ( 18)). One downside of the approach is that it includes multiple penalization terms. Therefore, extra additional hyperparameters need to be tuned. These hyperparameters have been set by trial and error. Typical behaviours when they are not set correctly is either all exponents going to zero, or à and B becoming proportional with a very high amplitude (see section 4.2.2 for the explanation).

Conclusion

This paper introduced a novel algorithm to find from a noisy dataset nonlinear input feature combinations which are optimal for physical modeling. They are obtained by solving a minimization problem based on the mutual information between some specific input combinations and the quantity to model. The algorithm has been tested on synthetic cases, showing very promising results. The last section has demonstrated that it can be modified results from the first LN network (q y ) (Bottom): results from the second LN network (q u ).

The results are close to the wanted expressions defined by equations ( 17) and [START_REF] Battiti | Using mutual information for selecting features in supervised neural net learning[END_REF].

to exploit dimensional information to form physically-relevant combinations of input feature. Hopefully, this type of approach could help provide better data-driven models for physics-related problems in the near future. As mentioned in the article, the framework could be modified to generate different forms of feature combinations that may be better adapted to some specific modeling problem, distinct from those considered in the paper. Note that section 2.2 has shown that some feature combinations that eliminate the symmetries of the model may provide significantly enhanced results over other reduced variables. As is, the algorithm does not account for possible symmetries in the data. This may be a future direction to explore to attempt to improve the algorithm. One identified downside of the approach developed here is that it involves several hyperparameters that need to be tuned adequately (but this is an issue for most of the existing data-driven techniques, unfortunately). The paper's results were relatively robust with respect to these hyperparameters.

But the synthetic cases considered were rather simple. It may be helpful to dedicate future studies to a more thorough analysis of the hyperparameter robustness in more complex cases. As the paper's primary goal was to introduce a new methodology for input feature design, exhaustive studies on the optimal choices of hyperparameters are left for future work, which may help to find empiric rules or automatic techniques to set them.

But arguably, the most interesting future work concerns the application of the FeDis technique to some of the open modeling problems mentioned in the introduction, such as RANS modeling for turbulent flows (open-source databases are already available to explore this path, see [START_REF] Mcconkey | A curated dataset for data-driven turbulence modelling[END_REF]). The mathematical and numerical framework introduced here is general and could be applied theoretically to any dataset. Nonetheless, it is not excluded that processing actual physical data involving advanced feature dependencies may raise implementation challenges that did not appear in this first study. For instance, since the size of the synthetic datasets from the present study was small, the question of the computational cost and convergence speed has not been investigated. Another interesting aspect is that actual non-synthetic datasets may involve "non-homogeneous" data stemming from multiple physical phenomena. Therefore, it may sometimes be hard to get a single model for the whole dataset, and designing multiple coexisting data-driven models may be better. Each of these models may have its own distinct set of relevant features. Therefore, exploring how the present approach may be coupled with clustering techniques may be an interesting question for future work.

Acknowledgment

This work is funded by ONERA as a part of the MODDA (MOdelization Data-Driven for Aerodynamics) project. state is dumped on disk and used for testing. This criterion is deactivated during the first fifteen epochs to avoid premature stopping.

Figure 1 :

 1 Figure1:(Left): Toy model f . The black square represents the range of values for (q 0 , q 1 ) used for training neural networks to estimate f (q 1 = 0 excluded).

Figure 2 :

 2 Figure 2: (Left): learning curves of N (naive strategy). (Right): learning curves of Ñ (model-informed strategy). The loss is the mean square error between the network output and the target value. Training is stopped based on a criterion defined in Appendix A.

Figure 3 :

 3 Figure 3: (Left): estimation of f produced by N (naive strategy). (Right): estimation of f produced by Ñ (model-informed strategy). The black square represents the range of values for (q 0 , q 1 ) used for training (everything outside the square shows extrapolation capabilities).

n

  i=1 a mi log(|q i |)

Figure 5 :

 5 Figure 5: Logarithmic Network (LN) structure. The hyper-parameter m is chosen depending on the number of input features chosen. In this article, given the strategy proposed, m will always be 1.

Figure 6 :

 6 Figure 6: Evolution of the learned exponents for each input variable.

Figure 7 :

 7 Figure 7: Evolution of the learned exponents for each input variable, second toy model (equation (14)). (Top): first run producing the first reduced variable. (Bottom): second run producing the second reduced variable, obtained by adding a penalization promoting new input feature combinations.

Figure 9 :

 9 Figure 9: Evolution of the learned exponents for the Law of the Wall test case. (Top):results from the first LN network (q y ) (Bottom): results from the second LN network (q u ). The results are close to the wanted expressions defined by equations (17) and[START_REF] Battiti | Using mutual information for selecting features in supervised neural net learning[END_REF].

  get the highest exponent m ← |a| > 0.1a max build a boolean mask for values lower than 10% of a max a ← a • m set all small values to zero ã ← nonzero(a) store all nonzero values a min ← min i (|ã i |) get the smallest nonzero exponent a ← a/a min normalize exponents end procedure

Description of the algorithm 4.2.1. Strategy

The algorithm is described in figure 8. It involves two LN networks, respectively producing a feature combination q y and q u of the form q y = y ∂u ∂y

q u = u ∂u ∂y

Note that the first weights A 0 and B 0 (the exponents of y and u respectively) are frozen to 1. Given the considered model (equations ( 17) and ( 18)), the expected outcome is A = (1, 0.5, -0.5, 0.5, 0, 0, 0, 0) and B = (1, -0.5, -0.5, 0.5, 0, 0, 0, 0). The rest of the algorithm is similar to the FeDis approach: the exponents are computed by maximizing a lower bound of the mutual information I(q y , q u ), following the same approach as that of algorithm 1. The specificity here comes from the penalization terms detailed in the next section.

Penalization of the loss

In addition to the classical L1 penalization of the exponents (used for all results in the paper), two specific penalization terms have been added for this algorithm. The first one aims to promote non-dimensional combinations of the features. The dimensional matrices M d,y and M d,u associated with q y and q u , respectively, are presented in table 1. They gather the physical dimension of each feature. The dimension has been set arbitrarily for the dummy variables ω i . The dimension of the LN combinations corresponding to A and B is then simply given by the vector-matrix products D y = AM d,y and D u = BM d,u (yielding a 3 × 1 vector). The extra penalization is then based on the L2-norm of these vectors:

with λ d a weighing coefficient for the penalization.

The algorithm has an easy (but unwanted) way to maximize the mutual information between q y and q u . If the subvectors à = (A 1 , . . . , A 7 ) and B = (B 1 , . . . , B 7 ) become aligned and their components becomes really large, it may lead to q y ≈ q α u , yielding a very high mutual information I(q y , q u ). This Appendix A. Regression task: section 2.2

The networks are made of four layers with 150 hidden units, eLU activation functions, followed by a last linear layer to produce the output. The loss function is the standard mean square error between output and target values. The network is trained using mini-batches of 100 samples. The optimizer is based on the Adam algorithm (learning rate of 10 -4 ). Training is stopped when the validation loss exceeds the lowest value encountered so far by more than 2% for more than thirty consecutive epochs. Then, the model's current The function T ψ from the algorithm is a neural network made of four layers with 200 hidden units, eLU activation functions, followed by a last linear layer to produce the output (scalar value). The optimizer is a Stochastic Gradient Algorithm (SGD) with a learning rate of 7×10 -2 . As mentioned in the article, an L1 penalization on the weights of the LN is added, this penalization is weighted in the total loss using a discount factor λ L1 = 5 × 10 Appendix C. FeDis algorithm: section 3.3

The function T ψ from the algorithm is a neural network made of four layers with 200 hidden units, eLU activation functions, followed by a last linear layer to produce the output (scalar value). The optimizer is a Stochastic Gradient Algorithm (SGD) with a learning rate of 5 × 10 -2 . As mentioned in the article, an L1 penalization on the weights of the LN is added, this penalization is weighted in the total loss using a factor λ L1 = 5 × 10 -3 . Training is performed using mini-batches of size 200.

For the second run of the algorithm, the extra penalization (alignment penalization, equation ( 16)) is weighted in the total loss using factor λ a = 3 × 10 -2 . All other hyperparameters are unchanged.

The losses evolution, not shown in the paper, is visible in figure C.11. The function T ψ from the algorithm is a neural network made of four layers with 200 hidden units, eLU activation functions, followed by a last linear layer to produce the output (scalar value). The optimizer is a Stochastic Gradient Algorithm (SGD) with a learning rate of 2 × 10 -2 . As mentioned in the article, an L1 penalization on the weights of each LN is added, this penalization is weighted in the total loss using a factor λ L1 = 5 × 10 -2 . Training is performed using mini-batches of size 400.

The first extra penalization (dimensionless promotion [START_REF] Bolón-Canedo | A review of feature selection methods on synthetic data[END_REF]) is weighted in the total loss using factor λ d = 5 × 10 -2 . The second extra penalization (alignment penalization, equation ( 16)) is weighted in the total loss using factor λ a = 1.

The loss evolution, not shown in the paper, is visible in figure D.12.


J F Semblat
email: semblat@lcpc.fr

P Dangla

M Kham

A M Duval

. The amplification of seismic motion is analysed in terms of level, occuring frequency and location. For both sites, the amplification factor is found to reach maximum values of 20 (weak motion). Site effects nevertheless have very different features concerning the frequency dependence and the location of maximum amplification. For the shallow deposit in Nice, the amplification factor is very small for low frequencies and fastly increases above 1.0 Hz. The irregular Caracas basin gives a much different frequency dependence with many different peaks at various frequencies. The model for Caracas deep alluvial basin also includes a part of the local topography such as the nearest mountain. One can estimate seismic site effects due to both velocity contrast (between the basin and the bedrock) and local topography of the site. Furthermore, the maximum amplification is located on the surface for Nice, whereas some strong amplification areas also appear inside the basin itself in the case of Caracas. One investigates the influence of this focusing effect on the motion vs depth dependence. This is of great interest for the analysis of seismic response of underground structures. The form and the depth of alluvial deposits are then found to have a great influence on the location of maximum amplification on the surface but also inside the deposit for deep irregular basins. It is essential for the analysis of the seismic response of both surface and underground structures.

Introduction

The analysis of seismic site effects considers amplification versus frequency curves showing the range of the spectrum leading to large motion amplification. Experimental measurements are generally performed along the surface with various methods : microtremor recordings, real earthquakes measurements [START_REF] Duval | Relation between curves obtained from microtremor and site effects observed after Caracas 1967 earthquake[END_REF]. Information on in-depth motion could sometimes be obtained thanks to specific measurement networks [START_REF] Kashima | Underground earthquake recording at Kushiro JMA observatory. 9th World Conference on Seismic Zonation[END_REF][START_REF] Lussou | Seismic design regulation codes: contribution of K-net data to site effect evaluation[END_REF]. Through numerical methods, one can also study the amplification process in various types of geological structures. It is for instance possible to consider the vibratory resonance of alluvial basin [START_REF] Bard | The two-dimensional resonance of sediment filled valleys[END_REF][START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF]. Otherwise, one can perform numerical analyses on site effects through explicit wave propagation models. In this paper, we try to study the influence of the basin geometry on site effects. Both surface and in-depth motion are especially considered to find out how they can be modified by some specific motion amplification for a typical basin geometry. The focusing effects are for instance taken into account to explain the possible increase of in-depth motion in some areas [START_REF] Sommerville | Seismic hazard evaluation. 12th World Conf. on Earthquake Eng[END_REF]. To perform such an analysis, seismic wave amplification is investigated in various types of alluvial basins considering the boundary element method.

Shallow and deep alluvial basins

For the analysis of in-depth motion amplification and focusing effect, we chose two alluvial basins with very different profiles : the first one is located in the centre of Nice (France) and is a wide flat basin (width 2 km, depth 60 m) [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF], the second one, located in Caracas (Venezuela), is a deep irregular valley surrounded by mountains (width 3.6 km, depth 300 m) [START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF]. Some experimental or numerical investigations were performed previously for both basins [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF][START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF][START_REF] Duval | Relation between curves obtained from microtremor and site effects observed after Caracas 1967 earthquake[END_REF]. We found that the amplitude versus frequency dependence is very different in each case. We will then try to analyse the variations of in-depth motion in both cases and to find out if focusing effects can actually influence motion amplification in a deep irregular basin.

Modelling site effects by the BEM

The numerical analysis of site effects for both types of basin was performed by the Boundary Element Method [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF][START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF][START_REF] Bonnet | Boundary integral equation methods for solids and fluids[END_REF][START_REF] Dangla | A plane strain soil-structure interaction model[END_REF]. The method is very powerful since it allows the modelling of seismic wave propagation for large geological structures without such drawbacks as numerical dispersion for some other methods [START_REF] Semblat | Efficiency of higher order finite elements for the analysis of seismic wave propagation[END_REF].

The numerical analysis was performed considering plane seismic waves of various types [START_REF] Dangla | A plane strain soil-structure interaction model[END_REF][START_REF] Betbeder-Matibet | In-depth attenuation of seismic ground motion[END_REF]. The shear wave velocities were chosen as follows: for Nice [1] C 1 =300m/s in the deposit and C 2 =1400m/s in the bedrock ; for Caracas [START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF] C 1 =450m/s and C 2 =2500m/s respectively. Fig. 1 gives the isovalues of the amplification factor for both sites. The first one (Nice, French Riviera) is shallow and its geometry is very regular. Site effects are found to be strong in the deepest part of the deposit (left) between 1 and 2 Hz and in the thinnest part (right) for frequencies above 2 Hz [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF]. For the second one (Caracas, Venezuela), there is a significant influence of the local topography (nearest mountains) as shown in Fig. 1 for 0.6 Hz [START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF]. The surface motion amplification has a very complex dependence on frequency. The irregular form of this basin as well as the large velocity contrast suggest that focusing effects could occur in the basin itself and influence the amplification process. In the following, we will estimate in-depth motion variations to determine if they can lead to deep amplification areas. This issue is very important for the design of earthquake resistant underground structures [START_REF] Semblat | Amplification and diffraction of seismic waves from underground structures[END_REF].

Occurrence of focusing effect

Focusing effect is related to particular geological structures that can focus seismic energy because of their geometrical and mechanical features. Some unexpected localised zones of damage were especially observed after Northridge earthquake [START_REF] Sommerville | Seismic hazard evaluation. 12th World Conf. on Earthquake Eng[END_REF]. To analyse potential focusing effect, we compare the seismic motion amplification at various depths for both alluvial basins (shallow regular ; deep irregular). Fig. 2 displays the amplification factor in the whole shallow basin (Nice) at various frequencies. It is given versus depth and distance (along the free surface). For the lowest frequency (1.0 Hz), there is only one amplification area on the free surface and in-depth motion decreases regularly. For frequencies values of 1.4 and 1.6 Hz, there are several amplification areas along the free surface in the left deepest part of the basin. No significant amplification is observed in the thinnest part (right) for those frequency values. For larger frequency values (2.0, 2.2 and 2.4 Hz), many different amplification areas are obtained along the free surface in the left part of the basin except for the last frequency value leading to low amplification in this part. In the right part, there is a strong increase of the amplification factor values for the three largest frequencies. Nevertheless, the seismic motion amplification is always decreasing with depth inside the basin. For all frequencies, there is a monotonic decrease of ground motion values from maximum surface motion values. Since the shallow basin is very flat, no focusing effects is observed but there is still a basin effect leading to amplification values much larger than those obtained from 1D analytical estimation considering the mechanical features of the deposit [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF]. For the shallow regular basin in Nice, site effects are then influenced by basin effects leading to seismic waves trapped in the deposit. However, there is no energy focusing effect due to the basin geometry and consequently no large in-depth amplification.

In the case of the deep irregular basin in Caracas [START_REF] Semblat | Seismic site effects in a deep alluvial basin: numerical analysis by the boundary element method[END_REF], amplification values versus depth and distance are given in Fig. 3. For both first frequencies (0.4 and 0.8 Hz), there are one or several (respectively) amplification areas along the free surface and in-depth motion decreases regularly down to the bedrock. For the second value (0.8 Hz), the deep deposit also appears more sensitive to some basin edge effects than the shallow basin. For the third frequency (1.2 Hz), we can suspect some little focusing effect since there is a very slow decrease of in-depth motion on the right part of the basin. At the bottom of the deepest part of the basin, there is a rather large value of seismic amplification. The focusing effects is much clearer for frequency 1.4 Hz : in the deepest part of the basin there is a strong increase of in-depth motion. It corresponds to an area of strong motion amplification located inside the alluvial deposit. For larger frequencies (1.8 and 2.0 Hz), there are several parts of the basin where in-depth motion increases. For some places, deep amplification can reach similar values to those obtained along the free surface. At 1.8 Hz, three main areas lead to in-depth motion increase and there are six of them at 2.0 Hz. These results (Fig. 3) show a strong influence of focusing effects on in-depth motion amplification. In the next section, we will discuss the dependence of seismic motion on depth by comparing in-depth motion curves for this deep site at various frequencies.

Fig. 2 : Amplification in the whole shallow deposit (Nice) at various frequencies.

Fig. 3 : Amplification in the whole deep deposit (Caracas) at various frequencies.

Influence on in-depth motion

In the case of the shallow basin, there is always a regular in-depth motion decrease in agreement with the classical rules (Fig. 2). For horizontally multilayered media, a simple analytical analysis leads to explicit decreasing laws for in-depth seismic motion [START_REF] Betbeder-Matibet | In-depth attenuation of seismic ground motion[END_REF]. For a two-dimensional shallow regular basin, results given in Fig. 2 follow the same trend than analytical results in the multilayered case.

To investigate the influence of focusing effects on in-depth motion for the deep irregular basin, several curves giving seismic motion versus depth are considered (Fig. 4). The variations of in-depth motion are very different. In some places, there could be a strong increase of in-depth motion due to focusing effect. As shown in Fig. 4, at 1.4 Hz, there is a maximum of the seismic motion inside the basin in its deepest part. For larger frequencies (1.8 and 2.0 Hz), the wavelength is shorter and in-depth motion local maxima appear in other parts of the deposit with different energy focusing processes. For the largest frequency (Fig. 4), there are even several different large motion areas along the maximum depth. In Fig. 4, the local maximum is shown to appear between 200 and 250m. For large frequency values, there are then several related focusing effects corresponding to the focus of seismic waves in shallower areas of the basin at shorter wavelengthes. The focusing effects are then influenced by the geometry of the basin as well as the depth/wavelength aspect ratio.

Fig. 4 : In-depth motion at various locations and frequencies for the deep deposit (Caracas).

Conclusion

The analysis of seismic site effects for two very different alluvial basins (shallow, deep) gives interesting results on potential energy focusing effects [START_REF] Sommerville | Seismic hazard evaluation. 12th World Conf. on Earthquake Eng[END_REF]. For the shallow regular basin considered in Nice, there is no focusing effect and larger amplification is obtained along the free surface. The influence of the basin geometry, vs wavelength, is only observed on the location of maximum amplification areas (deepest part for low frequencies and thinnest part for higher frequencies). For the deep irregular basin in Caracas, various amplification areas are observed inside the basin itself starting in the deepest part of the basin at some intermediate frequency. For larger frequencies (shorter wavelengths), different parts of the basin lead to large deep amplification. In-depth motion variations are consequently influenced by focusing effects. The shallow basin gives a classical decrease of in-depth motion whereas the deep basin can lead to some in-depth motion increases due to energy focusing effects. It is of great interest for the design of earthquake resistant underground structures [START_REF] Semblat | Amplification and diffraction of seismic waves from underground structures[END_REF] as well as the analysis of seismic hazard in urban areas [START_REF] Guéguen | From soil-structure to site-city interaction[END_REF].

Fig. 1 :

 1 Fig.1 : BEM modelling of site effects for shallow and deep alluvial deposits : amplification factor in the case of Nice (top) and Caracas (bottom)

Seismic Site Effects for Shallow and Deep Alluvial Basins:

In-Depth Motion and Focusing Effect