Patient Ntumba
email: patient.ntumba@inria.fr

Nikolaos Georgantas
email: nikolaos.georgantas@inria.fr

Vassilis Christophides
email: vassilis.christophides@ensea.fr

Scheduling of Continuous Operators for IoT edge Analytics with Time Constraints

Keywords: IoT (Edge), data stream, Fog, Cloud, continuous operator, scheduling, real time, queuing model

Data stream processing and analytics (DSPA) engines are used to extract in (near) real-time valuable information from multiple IoT data streams. Deploying DSPA applications at the IoT network edge through Edge/Fog architectures is currently one of the core challenges for reducing both network delays and network bandwidth usage to reach the Cloud. In this paper, we address the problem of scheduling continuous DSPA operators to Fog-Cloud nodes featuring both computational and network resources. We are paying particular attention to the dynamic workload of these nodes due to variability of IoT data stream rates and the sharing of nodes' resources by multiple DSPA applications. In this respect, we propose TSOO, a resource-aware and time-efficient heuristic algorithm that takes into account the limited Fog computational resources, the real-time response constraints of DSPA applications, as well as, congestion and delay issues on Fog-to-Cloud network resources. Via extensive simulation experiments, we show that TSOO approximates an optimal operators' placement with a low execution cost.

I. INTRODUCTION

The Internet of things (IoT) delivers at real-time massive amounts of data. To support (near) real-time control and automation applications such as smart transportation [START_REF] Tiwari | Reconfigurable streaming for the mobile edge[END_REF] or security [START_REF] Dangal | Towards industrial security through real-time analytics[END_REF], distributed IoT data streams need to be processed and analysed on the fly. Thus, data stream processing and analytics (DSPA) engines are used to extract in real-time valuable information from unbounded data streams via a series of continuous operators (i.e. aggregation, filter, join, etc.) [START_REF] Röger | A comprehensive survey on parallelization and elasticity in stream processing[END_REF].

Figure 1 highlights a DSPA application for country-wide traffic monitoring, as well as, city-wide traffic regulation. The latter enables control of traffic lights [START_REF] Tiwari | Reconfigurable streaming for the mobile edge[END_REF] e.g., to give priority to jammed traffic flows over non-jammed ones. We need to specify DSPA operators that analyse the data streams produced by the moving vehicles (e.g., vehicle identifier, GPS-location, driving speed). These operators are usually deployed at the Cloud where IoT data streams are transmitted e.g., by 5G based street antennas via several hops of the wide area network (WAN) [START_REF] Deng | Optimal workload allocation in fog-cloud computing toward balanced delay and power consumption[END_REF]. At high data stream rates, the WAN bandwidth to reach the Cloud can become a bottleneck as it may incur high network delays that may impair on the real-time constraint of DSPA applications. Furthermore, the WAN bandwidth is inherently dynamic due to the varying internet traffic conditions [START_REF] Jonathan | Wasp: wide-area adaptive stream processing[END_REF]. When the sum rate of transmissions of traffic data from vehicles are overwhelming the bandwidth capacity agreed between the Cloud provider and ISPs, the monetary cost charged for the WAN bandwidth usage can be much larger than the cost for computational resource usage [START_REF] Amazon | Aws whitepapers guides[END_REF].

Deploying DSPA operators closer to where IoT data streams are produced through Edge/Fog computing is currently one of the core challenges for reducing both network delays and bandwidth compared to DSPA applications deployed exclusively on the Cloud [START_REF] Bittencourt | The internet of things, fog and cloud continuum: Integration and challenges[END_REF]. The Edge/Fog computing introduces intermediate layers of computing, (storage), and networking between the IoT devices and the Cloud. The two terms are often used interchangeably. However, according to [START_REF] Iorga | Fog computing conceptual model[END_REF], we consider them as different layers: Edge computing is the network layer encompassing the IoT devices and their users, to provide, for example, local computing capability on a sensor (e.g. smart phone, Raspberry Pi, connected vehicle, drone, etc.). Fog computing residing between IoT devices network (Edge) and central Cloud, consists of Fog nodes which are either physical or virtual components (e.g. routers, small server, cloudlet, etc.). However, Edge/Fog nodes come with limited resources both in terms of computational power (e.g., CPU/GPU, RAM, etc.) and power supply (e.g., batteries, etc.). Thus, we need to schedule DSPA operators across Edge-Fog-Cloud resources by controlling both the processing latency and the network delays while keeping minimal the cost of using these resources [START_REF] Ali-Eldin | The hidden cost of the edge: A performance comparison of edge and cloud latencies[END_REF]. Fig. 1: Traffic management application adapted from [START_REF] Tiwari | Reconfigurable streaming for the mobile edge[END_REF] For instance, in our motivating application, traffic data rates may significantly vary per geographic area due to the high variability of the traffic density, as well as, of the mobility patterns. Clearly, workload picks require a higher usage of computational resources, energy and bandwidth and may incur higher network delays and processing latency [START_REF] Ali-Eldin | The hidden cost of the edge: A performance comparison of edge and cloud latencies[END_REF]. Therefore, a dynamic scheduling at run time is needed in order to decide how/when to schedule DSPA operators w.r.t. the evolution of the workload [START_REF] Hung | Dynamic workload allocation for edge computing[END_REF]. In our previous work [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF], we introduced a cost model of both network (bandwidth, delay) and computational resources (CPU/RAM) of network nodes. Using this model, we formulated the problem of scheduling DSPA operators between the Fog and the Cloud nodes to process data streams produced by mobile IoT devices at the Edge as a single objective optimization (SOO) problem of the combined resource cost subject to the constraints of the computational resource usage, cloud bandwidth usage and operator replicability constraint. As SOO is a NP-hard problem, we introduced a greedy algorithm that approximates an initial optimal placement of operators to replicate from the Cloud to the Fog. However, SOO did not address dynamic re-scheduling of continuous operators that also satisfies real-time response constraints.

In this paper, we extend SOO with real-time response constraints, that we formulate as a Time based Single Objective Optimization (TSOO) problem. We show that this problem is NP-hard and hence, the contributions of this paper are the following:

• we formalize TSOO to dynamically allocate both computational (i.e., CPU/RAM) and network resources (i.e., network delays and bandwidth); • we introduce a framework for monitoring the resource usage of network nodes as well as the workload of a DSPA application in order to trigger dynamic operator scheduling in the Edge-Fog-Cloud continuum; • we propose a heuristic algorithm called TSOO to approximate the best feasible scheduling of operators at run-time. We experimentally evaluate TSOO by using the simulation tool iFogSim [START_REF] Gupta | ifogsim: A toolkit for modeling and simulation of resource management techniques in the internet of things, edge and fog computing environments[END_REF] while considering heterogeneous computational and network resources across Edge-Fog-Cloud architecture and IoT data streams with varying rates. As baseline we consider TRCS [START_REF] Ntumba | Scheduling continuous operators for iot edge analytics[END_REF] that improves a pure cloud based processing by using as less as possible Fog resources to ensure Cloud bandwidth usage constraint and the real-time response constraint. We also compare TSOO with state of the art methods based on simulated annealing (SA) (TSOO-SA) [START_REF] Delahaye | Simulated annealing: From basics to applications[END_REF]. We use SA because of its flexibility and global optimization capabilities to approximate the global optimum of NP-hard problems. Furthermore, we use the heuristic algorithm proposed by Rizou et al. [START_REF] Rizou | Fulfilling end-to-end latency constraints in large-scale streaming environments[END_REF] that minimizes only the network resource usage cost while satisfying a real-time response constraint.

Our simulation results show that TSOO approximates well (with only 5% difference ratio) the optimal solution produced by TSOO-SA of using computational and network resources both statically and at run-time. However, TSOO-SA produces an operator placement which is not likely to satisfy real-time response constraints. Unlike TSOO-SA solutions that may end-up with 100% of constraint violations, TSOO solutions violates only up to 6.6% the real-time response constraint in the worst case of high rate data streams.

The remainder of this paper is organized as follows: Section II introduces our core model of DSPA applications and formulates the resource allocation problem. Section III extends our core model with real-time response constraints that DSPA applications have to satisfy. Section IV presents our heuristic algorithm TSOO for approximating optimal placements of DSPA operators. Section V details the experimental evaluation. Finally, Sections VI and VII present respectively the related work and the conclusions.

II. BACKGROUND

In this section we present the core model of Edge-Fog-Cloud architectures and DSPA applications and we formulate the resource allocation problem (for more details, readers are referred to [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF]).

A. Edge-Fog-Cloud Architecture

We abstract an Edge-Fog-Cloud architecture as a hierarchical WAN of resources defined by the set H={E, F, C} [START_REF] Li | Online workload allocation via fog-fog-cloud cooperation to reduce iot task service delay[END_REF]. The Edge (E) layer consists of M IoT devices E={E1, ..., EM } moving in N geographic areas Aj, j={1...N }, the Fog (F) layer consist of N Fog nodes F={F1, ..., FN } where each Fog node Fj provides nearby computational service to the geographic area Aj. One Cloud node C is considered at the top of the hierarchy.

Edge-Fog-Cloud nodes are characterized by the following parameters: (i) cmE i the available CPU/memory of an Edge node Ei and cmF j of a Fog node Fj; (ii) ndE i F j the network delay of an Edge node Ei to its nearest Fog node Fj; (iii) nbEF j the available network bandwidth capacity of all incoming links to the fog node Fj; (iv) ndF j C the network delay of the link from a Fog node Fj to the Cloud node C; (v) nbF j C the available network bandwidth capacity on the Fog node Fj to Cloud network link. The CPU/memory of the Cloud node C, cmC , is practically considered unlimited [START_REF] Li | Online workload allocation via fog-fog-cloud cooperation to reduce iot task service delay[END_REF].

We additionally consider the usage of CPU/memory cmuE i , cmuF j and cmuC respectively on a Edge node Ei, a Fog node Fj and the Cloud node C, as well as the usage of the network bandwidth nbuE i F j and nbuF j C on links respectively from an Edge node Ei to a Fog node Fj and from a Fog node Fj to the Cloud node C. Finally, let Sj be the sum of data streams arriving to a Fog nodes Fj and produced by mj(t) IoT devices moving at a time t in a geographic area Aj. Hence, M = N j=1 mj(t).

B. DSPA application model

We model a DSPA application as a directed acyclic graph (DAG) of operators, denoted by G , where the vertices represent the operators and the edges the data stream flowing between two operators [START_REF] Röger | A comprehensive survey on parallelization and elasticity in stream processing[END_REF]. G topology further includes the sources that produce the raw data streams Sj of rate |Sj| consumed by DSPA operators and sinks that capture the stream of the computed results. To cope with the infinite nature of data streams, we consider that continuous operators are executed in time windows ωx to process a finite set of data items dx arising within a time interval.

The operators are characterized by the following parameters: i) selx: the operator selectivity defined as the ratio between the input and output data rate of an operator Ox; ii) λx,y: the edge data rate defined as the rate of data stream flow between two operators; iii) cx: the computational cost of an operator [START_REF] Djemai | Mobility support for energy and qos aware iot services placement in the fog[END_REF] defined in terms of CPU (e.g., million instructions (MI) demand per byte of data) and/or RAM (e.g., 1Mb demand per 10Mb of data) for an operator Ox to process a unitary data dx; and iv) Dx: the operator data load defined as the aggregation of the input data items dx per time window ωx:

Dx= I i=1 ωx • λi,x.
Where I is the maximum number of upstream operators Ox producing data stream at rate λi,x towards the operator Ox. Finally, we denote by reqx = Dx•cx the computational resources CPU/memory required by an operator Ox given its data load Dx at a time t and its cost cx.

To replicate and migrate a part of G in different computational nodes ni = Ei|Fj|C of an Edge-Fog-Cloud architecture, we need to partition G in disjoint subgraphs,Gmigi, according to some workload criteria. The resulting graph to deploy is defined as follows:

G dep = ∀ n i ∈H Gmigi (1)
To specify a replication and migration point in G, we rely on the edge-cut algorithm which partitions G in two disjoint subgraphs [START_REF] Ntumba | Scheduling continuous operators for iot edge analytics[END_REF]. An edge-cut ecj contains the set of edges having one endpoint in each subgraph of the partition. Let |ecj| be the rate of an edge-cut ecj defined as the sum of edge data rates crossing this edge-cut.

C. Resource allocation problem

To capture the fact that computational resources allocated to a DSPA application at different layers of an Edge-Fog-Cloud architecture are limited, we weight the computational resource usage of different operators sharing a node by the inverse of the available computational resources at this node. Our model relies on the following intuitive assumptions [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF], [START_REF] Li | Online workload allocation via fog-fog-cloud cooperation to reduce iot task service delay[END_REF]: (i) the computational resources are practically unlimited in the Cloud, cmC → ∞, thanks to the on-demand resource scaling, and hence the weight in the cloud is practically zero, 1 cm C → 0; (ii) the computational resources of Fog and Edge nodes are limited as they can not be scaled on demand, and thus the weight in a Fog (or Edge) node is non-zero,

1 cm F j ∈]0, 1] (or, 1 cm E i ∈]0, 1
]). To simplify our model, we consider that mobile IoT devices at the Edge layer are not used for processing data streams [START_REF] Li | Online workload allocation via fog-fog-cloud cooperation to reduce iot task service delay[END_REF]. This will be address in the future work. Then, the focus of our work is to optimize the Fog computational resource usage cost:

cru = N j=1 cmuF j • 1 cmF j (2)
cmuF j is the sum of the RAM (or CPU/GPU) usage required (reqx) by each operator of a subgraph Gmigj, to replicate on Fj.

By abstracting the Edge-Fog-Cloud architecture as a hierarchical WAN of resources, two conflicting factors are involved [START_REF] Bittencourt | The internet of things, fog and cloud continuum: Integration and challenges[END_REF]: (i) the network bandwidth is increasing in the hierarchy as we go from the Edge to the Cloud and (ii) the network delay that also increases up as the Cloud is in two hops from the Edge and the Fog at one hop.

Unlike [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF], we consider the fact that the network delays and network bandwidths of each individual WAN links can be dynamic with regard to the network conditions [START_REF] Jonathan | Wasp: wide-area adaptive stream processing[END_REF]. Thus, we need to differentiate network links not only according to their delay [START_REF] Cardellini | Optimal operator placement for distributed stream processing applications[END_REF] but also by their available bandwidth [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF]. In this respect, the bandwidth usage of a network link is weighted by the inverse of its available bandwidth multiplied by the ratio between its network delay and the maximum network delay (ndmax) among all links in the hierarchical WAN of resources: ndmax = max i,j

{ndE i F j , ndF j C }.
As the Edge nodes are not used for processing data streams, each entire raw data stream Sj produced by an IoT area Aj at the Edge reaches the Fog anyway. Thus, the cost part concerning the Edge to Fog network links is fixed (set as c constant), we then focus on minimizing the Fog to Cloud network resource usage cost:

nru = c + N j=1 nbuF j C • 1 nbF j C • ndF j C ndmax (3)
The network bandwidth effectively used on all the Fog-to-Cloud network links is defined as follows: B = N j=1 nbuF j C

III. RESPONSE TIME AWARE RESOURCE ALLOCATION PROBLEM

In the core model of section II, we assumed that the operators of G can be deployed only between the Fog and Cloud. According to the criteria of minimizing cru and nru, the resulting G dep defined in Formula (1) becomes the disjoint partition in subgraphs Gmigj to deploy on the Fog nodes Fj and GmigC to deploy on the Cloud node C. In this respect, each Gmigj is delimited with GmigC by a replication and migration point, an edge-cut ecj. However, G dep should also take into account any real-time response constraints imposed to a DSPA application. Thus, prior to introduce the problem statement, we formalize the response time of DSPA application.

A. Response time model

We define the response time T as the worst end-to-end latency Lπij among all operator paths πij in G dep [START_REF] Cardellini | Optimal operator placement for distributed stream processing applications[END_REF], each data stream Sj is processed by nπ > 0 operator paths πij ∈ G dep with i={1 . . . nπ}:

T = max π ij ∈Gdep (Lπij) (4)
To calculate the end-to-end latency Lπij of an operator path πij, we consider the network delays of each network link traversed by this operator path along with the latency of each operator Ox ∈ πij for processing its data load Dx:

Lπij = max ∀E i ∈m j (t) (ndEiF j ) + exy ∈π ij nd M(Ox),M(Oy ) + Ox∈π ij lx (5)
1) Network link delays: Given that data streams are exclusively processed between the Fog and Cloud layers, Formula (5) includes only the worst network delay ndEiF j among all Edge to Fog network links. Then, M is the mapping function, that gives the resource node nj=Fj|C on which an operator Ox is (or can be) mapped to, and nd M(Ox),M(Oy ) is the network delay for transmitting data between two resource nodes on which the operators Ox and Oy are mapped to. nd M(Ox),M(Oy ) is negligible if the operators Ox and Oy are placed on the same resource node.

In general, the network delay includes the propagation delay of the network link that depends on the distance between the connected nodes, the processing and queuing delays of a packet at the intermediate routers and the transmission delay which depends on the available bandwidth on the network link [START_REF] Lohrmann | Elastic stream processing with latency guarantees[END_REF]. By assuming that the processing delay is a small constant as the router need to analyse only the data packet header to get its destination and send it to its destination, the network delay can be formulated as the sum of the propagation delay and the transmission delay [START_REF] Prasad | Bandwidth estimation: metrics, measurement techniques, and tools[END_REF]:

ndn j n k = pdn j n k + tdn j n k .
a) Propagation delay (pdn j n k ): is the time it takes to transmit bits of data between two connected nodes and it is independent on the data size [START_REF] Rizou | Fulfilling end-to-end latency constraints in large-scale streaming environments[END_REF]. However, it depends on the type of the network link medium, the distance between the connected endpoints and is limited by the speed of the light. To approximate the propagation delay between two nodes, we can use the Vivaldi algorithm [START_REF] Rausch | Synthesizing plausible infrastructure configurations for evaluating edge computing systems[END_REF].

b) Transmission delay (tdn j n k ): is the time for sending data on a network link from node nj to node n k . It depends on the size of data to transmit and the available network bandwidth. The latter is impacted by several factors including the number of active sessions, the transmission capacity of the link (nominal network bandwidth capacity), the medium access control, etc. Several techniques have been proposed to estimate the available network bandwidth of a network link [START_REF] Prasad | Bandwidth estimation: metrics, measurement techniques, and tools[END_REF]. In our work, to estimate the transmission delay of any data dxy of size |dxy| from the operator Ox mapped on the node nj to the operator Oy mapped on the node n k , we first evaluate the transmission delay td ′ n j n k of a data d of considerable size |d| from node nj to node n k , the actual transmission delay of data d is td ′

n j n k = td ′ n j n k -pd ′ n j n k .
Then, the transmission delay of any data dxy is calculated as follows [START_REF] Rizou | Fulfilling end-to-end latency constraints in large-scale streaming environments[END_REF]:

tdn j n k = td ′ n j n k • |dxy | |d|
2) Operator latency: The latency of an operator Ox depends on its current data load Dx, the type of operation it performs (e.g., filtering, projection, aggregation, etc.) and the available computational resource in terms of CPU (cpuj) of the hosting node nj = Fj|C. Thus, let µx be the rate at which an operator Ox can process its data load Dx on a node nj [START_REF] Renart | Distributed operator placement for iot data analytics across edge and cloud resources[END_REF]: µx = cpuj/reqx.

We assume that the resource nodes use a time sharing overbooking strategies to enable CPU allocation even if the CPU demand is greater than the total CPU capacity. Thus, if M IP Sj is the total CPU capacity of a resource node nj, then cpuj=min(M IP Sj, (M IP Sj/(q + qj))). Where, q and qj are respetively the number of current running processes and the operators to replicate on the node nj [START_REF] Rizou | Fulfilling end-to-end latency constraints in large-scale streaming environments[END_REF].

Given the infinite nature of data stream, let wtx be the waiting time of data elements in an operator queue. The service rate µx, the waiting time wtx and the number of data elements in an operator queue are random variables over a continuous time parameter. It is not so easy to calculate the operator latency lx. Thus, we model each operator as a queuing system with one server and following the first in first out policy [START_REF] Shortle | Fundamentals of queueing theory[END_REF]. Then, the operator latency lx is approximated as follows: E(lx) = E(wtx) + 1 µx . To approximate wtx, we need to consider the characteristic of each operator Ox in G, in particular, whether it relies on count based or time based windows.

Determining the characteristic of each operator enables to describe each operator queuing model with the Kendall's notation in the form of X/Y/1 where X is the data arrival rate distribution, Y the service rate distribution and 1 server [START_REF] Shortle | Fundamentals of queueing theory[END_REF]. Consequently, to formulate the associated waiting time E(wtx).

In this paper, we assume that the operators implement time based sliding window. This window consider the invariable temporal extend of the window called window time ωx, the progression temporal step, called sliding time βx where ωx > βx [START_REF] Jiang | Queueing analysis of relational operators for continuous data streams[END_REF]. The window contains the set of data that arrives within the last ωx time units, so that window data are processed every βx time units. The size of each window Dx is dynamic and dependents on the actual IoT data stream rate. The data arrival rate λx to an operator Ox may follow an exponential distribution [START_REF] Jiang | Queueing analysis of relational operators for continuous data streams[END_REF]. As the service rate depends on the size of the data to process, it also follows an exponential distribution. Hence, we model a time based sliding window operator as a M/M/1 queuing system where the approximated wtx is:

E(wtx) = λx µx•(µx-λx) .

B. Optimization Problem

Given the overall (Fog) computational resources usage cost, i.e., cru and the overall (Fog-to-Cloud) network resource usage cost, i.e., nru, the optimal placement M of Gdep, between the Fog and Cloud layers should minimize both cru and nru while satisfying the resource usage and the real-time constraints. As simultaneously minimizing nru and cru is conflicting, we consider a single objective optimization of the overall resource usage cost RU defined as the weighted sum of the normalized forms of cru and nru, respectively to CRU and to N RU by applying the min-max scaling technique:

minimize RU = wc • CRU + wn • N RU (6)
subject to cmuF j ≤ cmF j j = {1, . . . , N },

B ≤ Bmax (7)
T ≤ T max.

Where wc ≥ 0 and wn ≥ 0 are respectively the weights for computational and network resource usage cost, which enable to specify a usage preference between the two types of resources.

Formula (8) constrains the Fog to Cloud bandwidth usage by an upper threshold Bmax. Due to the limited computational resources in Fog, Formula (7) constrains the usage of these resources to guarantees that the CPU/memory usage of each Fog node Fj does not exceed the available CPU/memory. To ensure the real time constraint, Formula (9) imposes that the response time of a specific DSPA application should not exceed a threshold T max

For processing a data stream, each operator Ox ∈ G requires computational and network resources to process and transmit data stream with a certain latency (lx). We can show that the edges in G represent the precedence constraint between operators as we model G as a DAG of operators. Furthermore, by defining T as the maximum end-to-end latency Lπ ij among all the individual operator paths πij, we can show that the critical operator path is the path with the maximum end-to-end latency Lπ ij , with Lπ ij > T max. Hence, we need to deploy G (as G dep ) between the nodes ni = Fj|C to minimize the end-to-end latency of the critical operator path πij while ensuring fairness in resource usage [START_REF] Strusevich | Shop scheduling problems under precedence constraints[END_REF]. The latter is defined in terms of: (i) optimal trade-off between CRU and N RU , (ii) resource usage constraints (i.e. Formulas [START_REF] Bittencourt | The internet of things, fog and cloud continuum: Integration and challenges[END_REF] and ( 8)) and (iii) operator replicability constraint. This problem is an instance of the job shop scheduling (JSS) problem known to be NP-hard [START_REF] Strusevich | Shop scheduling problems under precedence constraints[END_REF]. Its complexity increases as we increase the number of resource nodes or the size of G.

Given that computational and network resources of H can be shared by several DSPA applications, the available resources may vary in time as long these applications can be deployed and removed on the fly. Moreover, the number and the rate of IoT data streams Sj may also vary according to the mobility patterns of IoT devices [START_REF] Hung | Dynamic workload allocation for edge computing[END_REF]. Under these conditions, an optimal placement M statically defined may not anymore be a feasible solution to our TSOO problem. For this reason, we need to dynamically reschedule an already deployed application graph G dep at run-time by identifying a new operator placement M that optimizes RU and satisfies the constraints.

IV. TIME BASED SINGLE OBJECTIVE OPTIMIZATION (TSOO)

To solve our TSOO problem, we need to search in G the edge-cut ecj of each data stream Sj, in order to build the graph G dep that includes each individual subgraph Gmigj to deploy on the Fog node Fj and GmigC to deploy on the Cloud node C. Each Gmigj should satisfy the constraint cmuF j ≤ cmF j and the operator replicability constraint. Moreover, the processed data streams Sj that will be transmitted on the Fog to Cloud network links should jointly satisfy the constraint B ≤ Bmax. Finally, the DSPA application response time T should satisfy the constraint T ≤ T max.

We observe that when selecting an edge-cut ecj as the replication and migration point of a data stream Sj, we can individually calculate the effect of this selection on RU . Thus, we assume that the computational and network resource usage costs can be split for each individual stream Sj: RU = N j=1 RUj, where RUj=CRUj+N RUj and RUj, CRUj and N RUj are respectively the contributions to RU , CRU , and N RU for processing Sj.

Our solution, called TSOO Algorithm, extends the SOO algorithm [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF]; we recall below the main functions that we use: (i) RUminCut() minimizes directly RU . To do so, it first considers Grep as the subgraph of G containing the set of operators that can be replicated on the Fog, so as to ensure the operator replicability constraint. Then for each data stream Sj, it identifies the subgraph Gmigj ∈ Grep delimited by the edge-cut ecj ∈ Grep that produces the minimum RUj while Gmigj satisfies the constraint cmuF j ≤ cmF j . (ii) DataMinCut() minimizes N RU . It processes like RU minCut(). However, it identifies the subgraph Gmigj ∈ Grep based on the edge-cut ecj ∈ Grep that produces the minimum N RUj.

A. Algorithm description

TSOO algorithm (cfr. Algorithm 1) starts by applying RU minCut() to generate a solution that attempts to minimize directly RU . If the output solution of RU minCut() satisfies the constraints B ≤ Bmax and T ≤ T max, then the achieved RU is optimal (line 6). Otherwise the problem may not have a solution or, if a solution exists, finding the optimal solution is NP-hard.

As a next step, TSOO algorithm applies a greedy search that produces local optimal solutions to approximate the global optimal solution in a reasonable amount of time. In this respect, we apply the function dataM inCut() to identify the solution that minimizes N RU and consequently B (lines 7-8). If the output solution of dataM inCut() does not satisfy the constraint B ≤ Bmax, the TSOO problem has no solution satisfying this constraint, unless we relax Bmax in order to accept this solution as the least bad one. If the constraint B ≤ Bmax is satisfied, TSOO further checks whether the constraint T ≤ T max is satisfied or not: (i) if the constraint T ≤ T max is satisfied, this means that dataM inCut() produces a good solution satisfying all the problem constraints; (ii) if the constraint T ≤ T max is not satisfied, starting from the solution of dataM inCut(), we search in the subgraphs Gmigj and GmigC the operators to move from the Fog to the Cloud (or the inverse) in order to reduce T until we satisfy the constraint (lines 10-13). This greedy search is performed with the functions operatorM oveback() and operatorM oveDown(), described in Section IV-B. If the greedy search produces a solution where the constraint is still not satisfied, the TSOO problem has no solution satisfying this constraint, unless we relax T max to accept this last solution as the least bad one.

If both constraints B ≤ Bmax and T ≤ T max are finally satisfied, TSOO applies a greedy search to minimize RU while keeping satisfied the problem constraints (lines 14-31).

Algorithm 1: TSOO Input: G, application graph Input: Grep, subgraph of replicable operators of G Input: Bmax, upper threshold for bandwidth usage Input: T max, upper threshold for response time Input: Sraw, set of raw data streams S j 1 X ← ∅, set of S j on which ∆RU is applied S j 2 M ← ∅, set of Gmig j , replicable subgraph per S j 3 RM ← ∅ set of edge-cuts ec j ∈ Grep per S j 4 RU ← 0, overall resource usage cost 5 B ← 0, Fog-to-Cloud network bandwidth usage 6 M ← RU minCut()

7 if B > Bmax then 8 M, RM, G dep ← dataM inCut(G) 9 if T > T max ∧ B ≤

B. Satisfy the Response time constraint

When moving an operator from the Cloud to the Fog, the operator latency will probably increase, assuming that the computational resources of a Fog node are smaller than the ones of the Cloud. At the same time, if other operators of the same operator path are already hosted on this Fog node, the latency of these operators will also increase, as the node resources will now be shared among more processes. Regarding the operators of the same path that remain in the Cloud after this move, we assume that the effect on their latency is negligible. On the other hand, the effect of this move on the network link delay between the Fog node and the Cloud depends on the size of the data produced by the moved operator in comparison to the size of the data that were transmitted from the Fog node to the Cloud before. In the opposite case, when moving an operator from the Fog to the Cloud, the operators latency on the same path will probably decrease, while the network link delay may increase or decrease.

In this respect, to satisfy the response time constraint we first apply the function operatorM oveback(), as it is more likely to reduce the response time T on an operator path. More specifically, we iteratively select the operator paths πij, where Lπij > T max, in decreasing order of their end-to-end latencies Lπij. For each selected πij, we start from the edge-cut ecj that delimits the two subgraphs Gmigj and GmigC , through which this operator path πij traverses, and we select the upstream replicated operator of ecj to be removed from the Fog node Fj. If this action improves the resulting end-to-end latency, we continue to remove the next upstream replicated operator, as long as the constraint B ≤ Bmax is satisfied. We stop applying operatorM oveback() if the constraint T ≤ T max is satisfied. However, if the constraint is finally not satisfied or if removing a replicated operator does not improve the resulting end-to-end latency, we next apply the function operatorM oveDown().

The operatorM oveDown() processes similarly. Rather than removing the replicated operators from the Fog, it replicates the not yet replicated operators on the Fog. Then, it stops if the constraint T ≤ T max is satisfied. Otherwise, the TSOO problem does not have a solution that satisfies the real-time response constraint. Due to space limitations, we omit the pseudo code of the two functions.

C. Run-time monitoring for rescheduling operator placement

DSPA engines such as Apache Flink [START_REF] Carbone | Apache flink: Stream and batch processing in a single engine[END_REF] run several DSPA applications. It consists of Job Manager (JM) and multiple Task Managers (TMs) distributed across the nodes. By assuming that data streams are processed on Fog and Cloud nodes, the JM can be deployed on Cloud as it provides practically unlimited resources and a TM can be deployed on each Fog and Cloud nodes [START_REF] Jonathan | Wasp: wide-area adaptive stream processing[END_REF]. The JM is responsible for planning the execution of operators by the mean of a scheduler and to assign them to the TMs. Each TM execute the assigned operators and it uses a local monitor to continuously monitors the operator metrics i.e. data arrival rate (λx) and service rate (µx), the computational resource usage (cmun i ) and the data stream size (nbuF j C ) to sent to the Cloud. Then it reports periodically these average monitored values to the JM for diagnosis. JM uses a global monitor to aggregate the reported values and to calculate: RU , B and T . and it decides whether M is still an acceptable solution under the following conditions:

RU ≤ RU max, ( 10
) kmax k=1 γ k = 0 ( 11
)
Where RUmax is an upper threshold of RU defined by the application owner and γ k a term that penalizes the violation of any constraint k of the TSOO problem under the following conditions:

γ k = max(0, val k -max k max k ) (12)
For each constraint k, val k can be the value of cmuF j , B or T and max k can be the value of cmF j , Bmax or Tmax.

If at least one of the conditions (10) and ( 11) is not satisfied, JM triggers the TSOO algorithm for rescheduling the operator placement. The current monitoring mechanism may incur overhead due to operator migration. Thus, existing mechanism for autonomous system such as Monitor, Analyze, Plan and Execute (MAPE) loop pattern can be used for this purpose [START_REF] Kephart | The vision of autonomic computing[END_REF].

V. PERFORMANCE EVALUATION

We use iFogSim to simulate an Edge-Fog-Cloud architecture as well as several DSPA applications. Besides the TSOO algorithm, our experimental framework also implements 3 competitor methods.

TSOO-SA: uses the simulated annealing (SA) meta-heuristic [START_REF] Delahaye | Simulated annealing: From basics to applications[END_REF] to solve our TSOO problem. Given that SA solves unconstrained optimization problems, we use a constraint relaxation method [START_REF] Delahaye | Simulated annealing: From basics to applications[END_REF]: we add to the cost function to minimize (i.e., RU ) a penalty ζ k for the violation of each individual constraint k of the TSOO problem:

f = wc • CRU + wn • N RU + kmax k=1 ζ k (13)
We set

ζ k = 1 if γ k > 0 (see Formula (12)), otherwise ζ k = 0.
SA starts from an initial solution or current solution fi and, based on an adequate perturbation method, generates a neighbour solution fj of the current one. SA accepts a neighbour solution to be the new current one on the basis of a probability controlled by a parameter called temperature τ . SA considers also: (i) the cooling rate α, which determines how quickly the temperature decreases; (ii) I, number of transitions for each value of τ , which determines after how many neighbour solutions generated, SA will decrease the temperature; and (iii) τ f inal , which determines the temperature at which SA will end. To set these parameters for TSOO-SA, we use the values proposed in [START_REF] Delahaye | Simulated annealing: From basics to applications[END_REF]: τ =300, τ f inal =10, α=0.9 and I=100 per each value of τ .

To generate a neighbour solution, we use the mutation technique Shift Move (SM) [START_REF] Hosseini Nasab | A simulated annealing heuristic for the facility location problem[END_REF]. In the current solution, we first represent as a key value table the mapping of each stream Sj and its replication and migration point ecj. Then, from this table, we randomly select a data stream Sj and we replace its current replication and migration point ecj by ec ′ j . ec ′ j is the next edge-cut of ecj in the source-to-sink direction of G. If ecj includes a sink of G, ec ′ j will take the value of the first edge-cut in G (source-to-sink direction). Then, we move to right for all the streams S k ̸ = Sj, the replication and migration points from one position.

TRCS (Time and Resource Constraint Satisfaction): extends the baseline approach RCS [START_REF] Ntumba | Scheduling continuous operators for iot edge analytics[END_REF]. TRCS enhances pure IoT Cloud analytics by dynamically placing the operators between the Fog and the Cloud in synergy with the evolution of the IoT data stream rates. More specifically, assuming an initial deployment of all the operators in the Cloud, TRCS minimally uses the Fog computational resources to satisfy the constraints T ≤ T max, B ≤ Bmax and B ≥ Bmin, where Bmin is a lower threshold of B used to avoid the oscillation of operator placement between the Fog and the Cloud.

Rizou et al [START_REF] Rizou | Fulfilling end-to-end latency constraints in large-scale streaming environments[END_REF]: propose an algorithm to distribute operators on a peer resource network. In this respect, the algorithm considers a latency (i.e., network delay) space where each host has a virtual position. It determines for each operator its optimal hosting node that minimises its network resource usage in the latency space, depending on the data stream rates of its neighbour operators. Starting from this solution, if the response time constraint is not satisfied, the algorithm moves each operator to the candidate host that increases the network usage as little as possible until satisfying the response time constraint.

A. Environment and system setup 1) Dynamic IoT data streams: We statically simulate the variability of the data stream rates arriving to the Fog nodes by selecting randomly (uniform distribution) 10 values of M IoT devices in the interval [5000, 50000], where each IoT device produces data at a rate of 4KB/s [START_REF] Xhafa | Evaluation of iot stream processing at edge computing layer for semantic data enrichment[END_REF]. Then, for each value of M, we set an interval [0, M ] in which we uniformly set mj(t) IoT devices per geographical area Aj, so that the sum of mj(t) be equal to M. In this respect, each Fog node Fj receives data stream Sj at rate |Sj| = mj(t)×4KB/s. We wish to have 15 results of T, B, RU, and CRU for each value of M and to plot the average of these results per value of M. Hence, we repeat the splitting of each value of M per geographical area 15 times. In this respect, the total data rate reaching the Fog follows a sequence of 10 uniformly distributed values of M ×4KB/s repeated 15 times. We feed this sequence to TSOO and the baselines.

2) Fixed Edge-Fog-Cloud architecture: The topology of the hierarchical Edge-Fog-Cloud architecture [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF] includes 1 Cloud node at the top of the hierarchy, 10 Fog nodes in the middle layer, and up to 50000 IoT devices at the Edge at the bottom. We use the tool Ether [START_REF] Rausch | Synthesizing plausible infrastructure configurations for evaluating edge computing systems[END_REF] to generate plausible network configurations of the Edge-Fog-Cloud architecture. The distribution of the resulting network configurations follows the one used in [START_REF] Jonathan | Wasp: wide-area adaptive stream processing[END_REF]. The propagation delay and network bandwidth for each Edge to Fog network links are respectively in the interval pdE i F j ∈ [10, 100] ms and nbEF j ∈ [START_REF] Hung | Dynamic workload allocation for edge computing[END_REF]50] Mbps and for each Fog to Cloud network links they are respectively in the interval pdF j C ∈ [100, 300] ms and nbF jC ∈ [100, 250] Mbps. For the computational resources, we simulate the Cloud node as an Amazon VM instance of type m6g.xlarge with 12GB of memory and 35900 MIPS of CPU [START_REF] Amazon | Aws whitepapers guides[END_REF]. At the Fog, we simulate ESXi virtual machines where the CPU capacities are set between 2400 and 8150 MIPS and the memory capacities are set between 1 and 4 GB [START_REF] Wöbker | Fogernetes: Deployment and management of fog computing applications[END_REF]. The MIPS evaluation comes from [START_REF]Zip LZMA Benchmark[END_REF].

3) DSPA application: We consider a model of the TLC application (New York City Taxi and Limousine Commission rides) [START_REF] Silva | Investigating edge vs. cloud computing tradeoffs for stream processing[END_REF]. The TLC application finds the busiest driver every two hours, where each vehicle emits at the end of a ride a data record containing driver identification, pick-up and drop-off times and locations. It comprises 5 operators and 1 sink (see Figure 2), where each operator uses a time based window and is modeled as a M/M/1 queue. Fig. 2: DSPA application used in the evaluation 4) Setting parameters: We set Bmax=125MB/s to avoid overloading the Fog to Cloud bandwidth usage, as in [START_REF] Ntumba | Efficient scheduling of streaming operators for iot edge analytics[END_REF] and for TRCS we set Bmin=30MB/s. Setting RUmax too high leads to high resource usage, which may lead to the system becoming unstable, while setting RUmax too low leads to inconsistent optimization of RU . Thus, we set RU max = 0.9. Given that the maximum propagation delay among all the network links is 300ms, we set T max=500ms. 2) Constraint satisfaction: As we need absolutely to have a solution to deploy for each value of M, we consider constraint relaxation when an algorithm does not find a solution that satisfies all the constraints. Therefore, we analyze whether an achieved RU satisfies or not our problem constraints. To do so, we divide by 15 the number of times, the value of B (or T ) fails to satisfy the related constraint. This gives the constraint violation rates, and shows how likely an algorithm can satisfy the constraints per each value of M.

Cloud bandwidth usage constraint (B ≤ Bmax). As depicted in Figure 4b, TSOO satisfies this constraint whatever the data stream rates for up to M =40000. However, TSOO does not satisfy this constraint at higher data stream rates (M ≥ 45000); it results in 100% of constraint violations. TSOO-SA fails to satisfy this constraint only for the highest data stream rate (M =50000 IoT devices), also with 100% of constraint violations. As for TRCS and Rizou et al., they do not satisfy this constraint starting from M ≥ 40000, with less than 100% of constraint violations.

Response time constraint (T ≤ T max). Figure 4c shows that TSOO satisfies this constraint whatever the data stream rates, except for the highest data rate (M = 50000), where it may fail to fulfill this constraint with 6.6% of constraint violations. TSOO-SA satisfies this constraint only at lower data stream rates (M ≤ 25000). Upon moderate or higher data stream rates, it may not satisfy this constraint with an increasing percentage of constraint violations (i.e. 6.67% to 100%) as the data stream rates increase. TRCS is likely to fail to satisfy this constraint with up to 93.3% of constraint violations upon M > 5000. This is due to the fact that TRCS prioritizes the satisfaction of the constraint Bmin ≤ B ≤ Bmax.

Finally, Rizou et al. has 0% of constraint violations for M < 40000, while it starts failing to satisfy this constraint for M ≥ 40000. However, the constraint violation rate is between 26.67% and 53.3%. This is attributed to the priority of the algorithm in satisfying the constraint T ≤ T max after identifying the minimum N RU .

3) Scalability analysis: We assess the scalability of TSOO against TSOO-SA and Rizou et al. in terms of execution cost, i.e., the time it takes for each algorithm to find the best operator placement, to rewrite the application graph based on this placement, and to deploy the resulting application graph. In this respect, we consider 5 different DSPA applications built based on the TLC application. Thus, App-1 has 5 operators (cf. Figure 2), App-2 has 6 operators, ..., and App-5 has 9 operators. We keep the same configuration of the Edge-Fog-Cloud architecture as presented in Section V-A2 and V-A4. We plot only the results for M=35000 IoT devices.

Figure 5 (A) shows that RU is practically constant for each one of the algorithms when identifying the operator placement of App-1, App-2 and App-3 with respectively 5, 6 and 7 operators. TSOO-SA produces the lowest RU among the three algorithms, and TSOO outperforms Rizou et al. As a matter of fact, the operator placement that each algorithm identifies for App-1 is the same for App-2 and App-3, at least regarding the resources that are taken into account in RU . On the other hand, TSOO-SA does not produce the lowest RU for the operator placement of App-4 and App-5, which have higher numbers of operators (8 and 9). This is due to the fact that the current operator placement, which is considered at run-time as the initial solution of TSOO-SA, may not be close to the optimal solution (as it is required [START_REF] Delahaye | Simulated annealing: From basics to applications[END_REF]). Thus, during the iterations for improving this solution, TSOO-SA fails to converge towards the optimal solution.

VI. RELATED WORK

Dynamic scheduling of DSPA applications is largely discussed in the context of resource rich execution environments(i.e. Cloud, cluster, etc.) than in the context of Edge/Fog computing due to the lack of solid established edge-oriented DSPA engines [START_REF] Cardellini | Run-time adaptation of data stream processing systems: The state of the art[END_REF], that why we use simulation tool to evaluate our proposed scheduling solution. In the context of resource-rich execution environments, Lohrmann et al. [START_REF] Lohrmann | Elastic stream processing with latency guarantees[END_REF] distribute DPSA application in cluster nodes with the objective to minimize resource usage cost and keep the response time under a certain threshold. Arkian et al [START_REF] Arkian | Model-based stream processing auto-scaling in geo-distributed environments[END_REF] propose a MAPE loop pattern for dynamic scheduling operators only at the Fog layer which providing high computational resources aiming to maximize the throughput. Heinze et al. [START_REF] Heinze | Latency-aware elastic scaling for distributed data stream processing systems[END_REF] focus on dynamic scheduling of DSPA applications in distributed homogeneous hosts. Then, their objective is to maximize the resource usage while keeping the response time under a certain threshold by using an operator migration technique. In our work, we are interested in continuously scheduling DSPA computations between heterogeneous Fog and Cloud nodes where the Fog nodes come with limited computational resources. Our objective is to satisfy both the response time and bandwidth usage constraints. In this respect, we are searching for an optimal trade-off between the Fog computational resource usage cost and the Fog to Cloud network resource usage cost.

Nowadays, machine learning (ML) techniques are used to recognize patterns from measured or profiled data to dynamically schedule operators in the Edge-Fog-Cloud continuum. For instance, reinforcement learning (RL) is used in [START_REF] Da Silva | Monte-carlo tree search and reinforcement learning for reconfiguring data stream processing on edge computing[END_REF] for dynamic scheduling that minimizes the DSPA application response time. The problem is modeled as a Markov Decision process (MDP) that explores two RL techniques (Q-learn and Monte-Carlo-Search-Tree). Furthermore, [START_REF] Russo | Reinforcement learning based policies for elastic stream processing on heterogeneous resources[END_REF] studies dynamic operator scheduling on heterogeneous infrastructure with the goal of minimizing the response time, resource usage cost and reconfiguration overhead. MDP model suffers from the dimensionality of the problem to model and hence the state space increases as the problem size increases that may bear high execution cost when comparing to a heuristic approach.

Meta-heuristic are largely used as ways of obtaining better solutions for NP-hard problems. For instance, SA is used in [START_REF] Fang | Sap: An iot application module placement strategy based on simulated annealing algorithm in edge-cloud computing[END_REF] to minimize the response time and energy consumption when placing application modules between the Cloud and the Edge. Particle swarm optimization is used by [START_REF] Djemai | A discrete particle swarm optimization approach for energy-efficient iot services placement over fog infrastructures[END_REF] to efficiently place application modules in the Edge-Fog-cloud continuum by minimizing energy consumption and response time while respecting the constraints of resource node capacities and module placement. The challenge in this context is to identify good setting parameters and a best scheme to generate neighboring solutions that may converge to an optimal solution. These approaches come with a high execution cost to be used for a dynamic operator placement.

VII. CONCLUSION

In this paper, we addressed the problem of dynamically scheduling operators between Fog and Cloud nodes, in response to the dynamic data stream rates produced at the Edge. We proposed TSOO, a resource-aware and time-efficient algorithm that takes into account the limited Fog computational resources, the real-time response constraints, and congestion and delay issues on Fog-to-Cloud network resources. Experimental results showed that TSOO is scalable and time efficient. It approximates the optimal solution by managing the trade-off between the usage costs of Fog computational resources and the Fog-to-Cloud network resources and satisfies the real-time constraint. As future work, we plan to extend TSOO to take as input the current operator placement from which a minimum set of changes will be applied in case of rescheduling rather than to run from scratch the whole algorithm. Furthermore, we plan to also include computational resources of the Edge layer.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Overall resource usage cost B. Evaluation results1)Resource usage cost: Figure3shows that TSOO approximates the optimal RU with an approximation error of up to 5.34%, when comparing to TSOO-SA. Furthermore, TSOO outperforms Rizou et al. and TRCS, whatever the number of IoT devices (and the produced data stream rates), with a difference ratio respectively up to 11.32%

Fig. 5 :

 5 Fig. 5: Resource usage cost (A) and Execution cost (B). In Figure 5(B), the execution cost is increasing with the number of operators per DSPA application. TSOO-SA has the highest execution cost, with a difference ratio of up to 729% when comparing to TSOO. TSOO's execution cost is very close to the one of Rizou et al. with a difference ratio of up to 0.12% (TSOO being slightly more costly).

  Bmax then 10 Set maxπ, sorted set of paths π ij where Lπ ij > T maxGet S j , e j k , Gmig j k corresponding to ∆ RU Rewrite G dep to include all Gmign j (or Gmig jp ) ∈ M 31 Send each Gmign j ∈ M to corresponding resource node n j

	17	Pull ∆RU on top of ∆RUset
	18	while ∆RU < 0 do
	19
	20	ec j ← RM [j]
	21	Calculate T when considering ec j
	22	B

11

operatorM oveBack(M, RM, G dep , maxπ)

12 if T > T max then 13 operatorM oveDown(M, RM, G dep , maxπ) 14 if B ≤ Bmax ∧ T ≤ T max then 15 ∆RUset ← edgeCutM ove(G, RM, M ) 16 Sort ∆RUset in increasing order ′ ← B -|ec j | + |ec j k | 23 if B ′ ≤ Bmax ∧T ≤ T max ∧ X.has(S j ) = False then 24 RM [j] ← e jp 25 M [j] ← Gmig j k 26 B ← B -|ec j | + |ec j k | 27 RU ← RU + ∆ RU 28 X.add(S j )

29

Pull ∆RU on top of ∆RUset 30


Audrey De Jong

Samir Jaber

Niall Ferguson

come


Hans-Beat Bürgi
email: hans-beat.buergi@unibe.ch

Alessandro Genoni
email: alessandro.genoni@univ-lorraine.fr

Remarks on X-ray Constrained/Restrained Wavefunction fitting Authors

Some comments concerning the meaning and significance of X-

Introduction

In 1998 Jayatilaka (1998) presented the "X-ray Constrained Wavefunction" (XCW) fitting approach, a method that combines quantum chemistry and X-ray diffraction to produce what the author called an "experimental wavefunction". The intention is to find a wavefunction and the corresponding electron density i) that is as close as possible to the crystal density , i.e. explains the observed structure factor amplitudes to within a weighted mean square deviation between experiment and model and ii) that simultaneously minimizes the quantum mechanical energy of the investigated system.

At present, XCW fitting (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;Grimwood & Jayatilaka, 2001;Bytheway, Grimwood & Jayatilaka, 2002;Bytheway, Grimwood, Figgis et al., 2002;Grimwood et al., 2003) is the last in a series of steps to model . The procedure usually begins with Independent Atom Modelling (IAM, with tabulated spherical atomic densities and form factors). It may be followed by interpreting the diffraction data with more refined descriptions of the atomic and bonding electron densities, such as Multipole Models (MMs) (Stewart, 1976;Hansen & Coppens, 1978) or quantum mechanical Hirshfeld Atom Refinements (HARs) [START_REF] Jayatilaka | [END_REF]Capelli et al., 2014;Kleemiss et al., 2021). They are usually based on Hartree-Fock (HF) or Density Functional Theory (DFT) calculations, with or without consideration of the environment of the structural fragment that is treated quantum mechanically [START_REF] Jayatilaka | [END_REF][START_REF] Wieduwilt | [END_REF]Ruth et al., 2022). HAR at correlated levels has also been attempted (Wieduwilt et al., 2020). At the last stage, the XCW procedure tries to extract information not included in the preceding steps, typically the polarization due to the crystal environment (Ernst et al., 2020) as well as correlation (Genoni et al., 2017) and relativistic effects (Bučinský et al., 2016).

The XCW procedure thus attempts to extract the last little bit of information from a diffraction experiment. In other words, it tries to improve an already quite detailed model of the electron density, whose R-factor is typically as low as 1-3%. The likelihood that this information is of similar magnitude as the systematic and random errors in the diffraction data and inadequacies of the quantum chemical model is thus not negligible and creates some unresolved problems concerning the meaning of the results obtained from XCW fitting. Some of these problems, not all of which have been addressed explicitly, are discussed below.

A summary of the X-ray constrained/restrained wavefunction approach

As mentioned in the introduction, the XCW fitting approach aims to determine a wavefunction that minimises the quantum mechanical energy of the considered system and simultaneously maximizes the agreement between experimental and calculated structure factor amplitudes (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;Grimwood & Jayatilaka, 2001;Bytheway, Grimwood & Jayatilaka, 2002;Bytheway, Grimwood, Figgis et al., 2002;Grimwood et al., 2003;[START_REF] Jayatilaka | Modern Charge-Density Analysis[END_REF]. The problem is mathematically equivalent to finding the wavefunction that minimizes the energy under the constraint , where is the squared goodness-of-fit given by the following expression:

is the r-th structure factor amplitude calculated from a model electron density ; is the r-th observed structure factor amplitude with as its estimated standard uncertainty; η is a scale factor that puts the calculated structure factor amplitudes on the same scale as the experimental ones; is the number of experimental observations, and the number of parameters. Furthermore, as indicated by Jayatilaka and collaborators in the seminal papers on the technique, in principle should be set equal to 1.0 (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;Grimwood & Jayatilaka, 2001;Bytheway, Grimwood & Jayatilaka, 2002;Bytheway, Grimwood, Figgis et al., 2002;Grimwood et al., 2003). This implies that the objective of XCW fitting is the determination of an "experimental wavefunction" that reproduceson average -the experimental structure factor amplitudes to within one standard deviation of the experimental measurements.

To minimize the energy of the system under the constraint expressed by relation (1), the following functional has been introduced:

The first term on the right-hand side is the quantum mechanical energy of the system (i.e.,

, with as the Hamiltonian operator), and is a Lagrange multiplier that has the dimension of an energy. To achieve the aim of XCW fitting, functional J must be made stationary with respect to both the wavefunction and the Lagrange multiplier .

However, in the current implementation of the XCW technique, the second stationary condition (i.e., the one with respect to ) is never considered and, for this reason, the functional that is really minimized is the following one: o This implies that the present XCW method simply attempts to determine the wavefunction that minimises the energy of the system under examination and reproduces as much as possible the measured X-ray diffraction data. thus loses its original meaning of a Lagrange multiplier; instead, it can be considered as the weight of the experimental data in the computations and simply becomes an external parameter that is gradually changed during the XCW procedure (see below). The term in equation ( 4) is a restraint in the crystallographic sense and the method should better be called "X-ray restrained wavefunction (XRW)" approach [START_REF] Jayatilaka | Modern Charge-Density Analysis[END_REF]Grabowsky, Genoni & Bürgi, 2017;Ernst, Genoni & Macchi, 2020;[START_REF] Macetti | [END_REF] or X-ray regularization procedure (Davidson et al., 2022a(Davidson et al., , 2022b)). For the sake of completeness, it is mentioned here that one of the present authors has very recently reformulated the Jayatilaka technique by explicitly considering the stationary condition of functional (3) with respect to (Genoni, 2022). Only in that case does regain its meaning as a Lagrange multiplier, and the approach can then be called "X-ray constrained wavefunction (XCW) method".

In this work, we focus on the original and current version of the Jayatilaka strategy, but refer to it as "X-ray restrained wavefunction (XRW) technique". urthermore, although multi-determinant versions of this method have already been developed (Genoni, 2017;Casati et al., 2017;Genoni, Franchini et al., 2018;Genoni, Macetti, Franchini et al., 2019), for simplicity and without losing generality, we limit our discussion to the original and most used version of the Jayatilaka strategy, namely the case in which the wavefunction to be determined has the analytical form of a single Slater determinant.

For the single Slater determinant wavefunction ansatz, it is possible to show (Jayatilaka, 1998;Jayatilaka & Grimwood, 2001;Grimwood et al., 2003) that finding the wavefunction which minimizes functional (4) corresponds to determining the molecular orbitals that satisfy the following modified Hartree-Fock equations (for a 2N-electron closed-shell system):

is the usual Fock operator of quantum chemistry, while the XRW operator is

Re Im with and . and in equation ( 7) are the real and imaginary parts, respectively, of the one-electron scattering operator where is the number of symmetry-unique positions in the unit-cell, are the rototranslations associated with the unit-cell symmetry operations ( and being a rotation matrix and a translation vector, respectively), is the reciprocal lattice matrix, and is the vector of Miller indices associated with the r-th reflection.

Note that the term can be seen as a perturbation with weight to the traditional Fock operator of quantum chemistry. However, depends on the experimentally determined structure factor amplitudes. Given that experimental measurements are subject to random measurement errors, is subject to random error too and will differ somewhat for different X-ray experiments performed on the same system under the same experimental conditions. In other words, not only accounts for shortcomings of the model of the electron density but is also characterized by a random component.

The interpretation of as a perturbation operator can be traced back to a pioneering idea by Weiss from the 1960s [START_REF] Weiss | X-ray determination of electron distributions[END_REF]. He envisaged the possibility of correcting the deficiencies of the Hartree-Fock model by recovering effects of electron correlation from the experimental X-ray diffraction data. More recently, the relation between Weiss' original ideas and the XRW approach has been highlighted by [START_REF] Macchi | Quantum Crystallography: Expectations vs Reality[END_REF] in a review on origins, state of the art and perspectives of quantum crystallography.

Expressing the molecular orbitals in equations ( 5) and ( 7) in terms of a set of basis functions (as it is usually done in quantum chemistry), the modified Hartree-Fock equations expressed by relation ( 5) straightforwardly transform to modified Roothaan-Hall matrix equations:

with as the typical Fock matrix of quantum chemistry calculations, C as the matrix whose columns contain the coefficients that expand the molecular orbitals in the chosen basis, as the diagonal matrix of the orbital energies, and as the matrix associated with the operator , which obviously has the same dimensions as the Fock matrix and is weighted by the external multiplier

In the current version of the XRW approach, equations (10) are solved self-consistently for different values of by following a "trial & error" procedure. This means that the molecular orbitals obtained through the computation with at the i-th step are used as guess for the calculation with . The procedure is usually iterated, but without a clear and definitive criterion to stop. Several different halting criteria have been proposed over the years. For a review see Davidson et al. (2022a); for two recent proposals see Davidson et al. (2022b). For reasons discussed in section 3, none of these seems convincing. The consequence is that although the XRW process should be halted at o , the final value of o often remains greater than the desired value or, under some circumstances, the calculations may go on beyond the desired limit to a final o value that is lower than 1.0.

The present version of the XRW procedure assumes that the nuclear positions and the atomic displacement parameters are known (to within their uncertainties). As also mentioned in the Introduction, these parameters may be obtained with the help of HAR (Capelli et al., 2014) since structural parameters resulting from Hirshfeld atom refinements of X-ray diffraction data are almost as accurate as those obtained through neutron diffraction experiments (Woinska et al., 2016). HAR also provides a starting point for the XRW procedure in the form of the final HAR wavefunction and, consequently, of the final electron density , which hereafter will be called and , respectively, with a corresponding squared goodness-of-fit . Note, however, that, given accurate atomic positions and ADPs, the quantities , and might also be obtained from a level of theory and a basis set different from those used for the Hirshfeld atom refinement.

No problem so far.

The subsequent XRW procedure modifies (and ) to (and ) until one obtains the best fitting of the experimental crystal density or, equivalently, the best fitting of its Fourier transforms . In other words, the XRW method tries to extract information from the measured structure factor amplitudes that is not already contained in , such as polarization, electron correlation and relativistic effects.

At least in principle, an XRW fitting can adjust all molecular orbitals coefficients. However, for a 2N-electron closed shell system described by a single Slater determinant, there are independent parameters to be determined (with as the number of basis functions used in the calculation). This number is usually about an order of magnitude larger than the number of observations . It is therefore important to choose the value or, alternatively, the target value for GoF 2 in such a way that the procedure produces the best approximation to , with as little contaminated as possible by imponderables of the diffraction experiment or the quantum mechanical method used, i.e. systematic and random error. 1 The problem of choosing has been termed the halting problem. For reasons 1 According to the IUCr definition, a systematic error is a contribution to the difference between an estimate and the true value of a quantity due to deficiencies of the model.

(https://www.iucr.org/resources/commissions/crystallographic-nomenclature/statdes/terms.html). In the case of structure factors, such contributions may arise from two sources. On the one hand, from an inadequate model of the diffraction experiment needed to relate the raw intensity observations to what is commonly called the observed structure factor amplitude, The symbol is somewhat misleading as it does not represent a direct experimental estimate, but rather a derived quantity that depends on assumptions and on estimates or measurements given in section 3, it is our opinion that no "tried-and-tested", convincing solution of this problem has been presented so far.

Some difficulties

To gain some more insights into the XCW fitting process, a few cases with specific assumptions about and are analyzed. Some of these assumptions are unrealistic but provide a reference for gaining insights into the fitting procedure applied to more realistic cases.

A) First of all, let us assume the unrealistic case with:

-, where is the random error in ;

-;

are trustworthy and normally distributed estimates of the random errors of ;

In this case, 2 and it is trustworthy. Therefore, no further fitting is indicated.

B) Making the same assumptions as in case (A) except that the uncertainties unknowingly underestimate the random errors of .

In this situation, 2 and further fitting is indicated. However, any such fitting would lead to wavefunctions and thus to structure factor amplitudes which include contributions from both and . In other words: for any non-zero value of λ, the XRW procedure merely fits experimental noise! For different experiments of the same substance under the same conditions, the random errors will be distributed differently and, consequently, the fit will look different. Underestimating is not uncommon in actual experiments. Note that this conclusion would not hold if, under the given assumptions, . However, this seems unlikely and has certainly not been shown to be the case.

of factors affecting the measurement process (background, X-ray absorption, X-ray beam polarization, extinction, etc.). If such factors are insufficiently accounted for, the derived quantity will be unreliable to an unknown degree. On the other hand, an inadequate model of the electron density distribution, e.g an insufficiently flexible basis used in the quantum mechanical calculation, can also cause systematic differences between and C) Another artificial case assumes that -(i.e., no random errors are included); -for all reflections;

is based on a quantum mechanical model with a basis set that is insufficiently flexible to reproduce .

This situation corresponds to the case of theoretically generated structure factor amplitudes and the fitting leads to and , where is the part of that cannot be modelled with . Given the shortcomings of the basis set, 0; 2 as well, but its expectation value is unknown because it depends on the unknown limitations of the model used.

2 is expected to decrease with increasing the number of basis functions (and thus the number of parameters that can be optimized).

For this kind of problem, a way to find an optimal has been proposed, albeit based on electron densities rather than on their Fourier transforms, the structure factors (Tozer et al., 1996). Equations analogous to equations ( 10) are solved by combining with three different values of λ and extracting an optimal λ from the three results. One of the procedures proposed in Davidson et al. (2022b), although inspired by Tozer et al. (1996) is different. It updates at each fitting step until a halting point is reached. It has not been shown that the two procedures are equivalent, i.e. lead to the same result. In any case, the important point is that the choice of affects the 2 independently of and in an unpredictable way.

D) For a more realistic case of XRW fitting, it is assumed that:

-; -;

is based on a quantum mechanical model with a basis set that is likely to be insufficiently flexible to reproduce ;

-the are trustworthy and normally distributed estimates of the random errors of .

In this case, the fitting leads to and , in analogy to case (C). However, in contrast to that situation, the modifications of in correct not only for inadequacies of , but also absorb into the model parts of , i.e. some of the experimental noise, similarly to case (B). The two contributions are of unknown magnitudes and thus has an undefined uncertainty on two accounts. Conversely accounts for those parts of and that could not be absorbed into . It is thus unclear again to know the expectation value of GoF 2 and .

E) The most realistic case in XRW fitting makes the same assumptions as in (D) except that:

-, where represents a systematic error of due to inadequacies of the processing of the primary data, e.g. the frames from the diffraction experiment;

-the are not trustworthy estimates of the random errors of because they are too small, too large or show an unusual (non-normal) distribution.

In this situation, the fitting leads to and , in analogy to cases (C) and (D). However, unlike (C) and (D), the calculated correct not only for inadequacies of , but also absorbs into the model parts of , i.e. some of the systematic and random errors. The three contributions are again of unknown magnitudes and thus has undefined meaning and uncertainty. Conversely accounts for those parts of , and that could not be absorbed into . If, in addition, the quality of the is unknown, estimating the expectation values of GoF 2 and is further aggravated.

In a recent paper, Davidson et al. (Davidson et al., 2022a) consider the XRW fitting approach as a regularization problem, i.e. the process of adding information in order to solve an ill-posed problem or to prevent overfitting [here we use the term 'regularization' as defined in https://en.wikipedia.org/wiki/Regularization_(mathematics)]. In the present case there are two ways of looking at this problem: either as adding experimental information to a quantum chemical model that is insufficient to determine on its own, or as adding theoretical quantum mechanical information to a model of the experiment that is usually strongly underdetermined because the number of coefficients in the wavefunction exceeds the number of observed structure factor amplitudes by far. In both cases the question arises as to how much extra information is required to get an optimal result, i.e. what is the optimal magnitude of λ and how should an optimal mix of calculated and experimental information be defined anyway? This is an alternative formulation of what has been called the halting problem (Davidson et al., 2022a;Davidson et al., 2022b). If one can reasonably assume that the model of equation ( 4) is sufficiently flexible to determine and that the are truthful estimates of the experimental uncertainties, one might choose a value of λ that results

in GoF 2 = 1. However, these conditions are rarely fulfilled implying that the choice of an appropriate final value for GoF 2 remains an open problem in the XRW scheme.

In both ways of looking at the problem, systematic and random errors also introduce uncertainty in the wave function, , not only in the choice of λ or of the final GoF 2 value.

Since the goal of the exercise is to find an "experimental" X-ray restrained wavefunction and thus an experimental electron density function , an additional question arises: what are the uncertainties in or ? Which parts of are reliable and which are not? The present XRW fitting algorithm provides no information on this question. In contrast, leastsquares modelling of the crystal electron density through more traditional methods (e.g., multipole models) determines a set of parameters, their standard uncertainties, and mutual correlations. There are statistical criteria for judging the significance of introducing additional parameters into the model, such as Hamilton's R-factor ratio test (Hamilton, 1965). Such criteria help preventing overfitting by the model. The uncertainty of can be estimated through propagation of error or analogous procedures. None of these options are available in the present XRW algorithm. Note that the problem raised in this paragraph would also affect the reformulation of the X-ray constrained wavefunction approach recently envisaged by one of the present authors (Genoni, 2022), although, in that case, if convergence is achieved, one can obtain the desired value Δ of statistical agreement between observed and calculated structure factor amplitudes.

In summary, there are several difficulties with the outcome of an XRW fitting and the meaning of an X-ray restrained wavefunction:

i) an X-ray restrained wavefunction models random error to an unknown degree; while the more traditional least squares optimization determines expectation values of parameters together with their uncertainties. XRW fitting with its quantum mechanical minimization condition determines parameters too, but not necessarily their expectation values and without uncertainty.

ii) an X-ray restrained wavefunction models systematic errors in the structure factors or inadequacies of the atomic basis set to an unknown degree (as does a least squares procedure), thereby obscuring the statistical meaning of GoF 2 .

iii) no procedure is available at present to quantify, at least approximately, the relative and absolute magnitudes of random and systematic errors, and of basis-set inadequacies in the optimized density .

iv) another non-negligible drawback associated with both the XRW and XCW approaches has been alluded to before: it is the problem of correctly defining GoF 2 and its expectation value. Even among developers of the XRW/XCW approach, it is still unclear today which value to assign to in equation ( 2). In some circumstances, it is simply set equal to the number of adjustable parameters used in the "wavefunction refinement" (such as when it is manually adjusted in the XRW calculations); in other cases, for example when the XRW computation follows an Hirshfeld atom refinement, includes the number of atomic positions and ADPs determined through the preliminary structural refinement; finally, in one of the seminal papers of the X-ray restrained wavefunction approach (Grimwood et al., 2003), it is pointed out that should vary as a function of , approaching the number of molecular orbital coefficients when becomes large (which in most practical applications would lead to a negative value in the denominator of the first factor on the right-hand side of equation ( 2)). This uncertainty in the definition of aggravates establishing the final and desired value of GoF 2 further and, consequently, determining the correct value of .

Therefore, due to above-mentioned problems, an X-ray restrained wavefunction loses -to an unknown degree -its intended meaning as a representation of , because it cannot be decided whether the XRW density , i.e. it cannot be decided whether or not the final electron density distribution resulting from the XRW process corresponds to the best possible representation of the real electron density . For the same reasons the halting problem is ill defined.

Possible remedies

What can be done to alleviate these uncertainties?

-The data: utmost care is needed to get minimally contaminated by systematic errors.

It should be clearly separated from diffuse scattering and very carefully corrected for the properties of the sample (absorption, extinction, etc.), the detector (non-linear response, oblique incidence, etc.) and other imponderables of the experiment. Such corrections are ideally based on physical rather than empirical models of the factors to be accounted for.

Good data, ideally coming from more than one sample, tend to minimize effects of some systematic errors, . Data should be tested for consistency between and within samples.

-The should be trustworthy and their distribution tested for consistency between and within samples. Experiments that produce small, trustworthy minimize effects of .

Good help to make GoF 2 quantities more meaningful, of course provided that the number of extra (effective) parameters coming from the XRW fitting can be estimated. So far this is not done. should be tested for outliers with defined criteria before and after XRW fitting, e.g. with normal probability plots.

-Separate analysis of experimental data should be done for different choices of data: splitting data into parts (cross validation; Krause et al., 2016), looking at different samples separately, and combining all data. It is expected that such a procedure leads to different because the influence of will be different. Comparisons and analysis of different for the same compound could be a first step towards characterizing experimental uncertainty.

-To gain a fuller picture of the influence of experimental uncertainty, including the part coming from the atomic positional and displacement parameters, the following iterative procedure might be used in validation experiments: the initial XRW electron density can serve to define Hirshfeld atoms which are used in turn for refining atomic positions and displacement parameters. This might be followed by another XRW fit, etc., until convergence is reached, following the spirit of the XWR (X-ray wavefunction refinement) approach already envisaged and proposed by Grabowsky and collaborators (Grabowsky et al., 2012;Woińska et al., 2017).

-Last but not least the results should be tested by analysing the same data starting from different models (e.g. HF, DFT and different basis sets). If the different XRW fittings give comparable results for , the analysis can be considered robust.

-One could also think of comparing results from alternative expressions for GoF 2 , e.g. one that includes the phases of the structure factors:

or some function of their Fourier transform, i.e. some function of the difference density.

Since the XRW fitting comes into play at the very end of the modelling, it can reasonably be assumed that the phases of and are the same (in the centrosymmetric case) or very close to each other (in the non-centrosymmetric case). This definition would be more akin to a restraint based on the difference electron density, provided one includes a complete set of structure factors up to the resolution limit of the experimental data and GoF 2 is divided by the volume V of the unit cell.

-For many of the open problems, a better understanding of the XRW procedure could be obtained from a well-crafted study with known, predefined , , , and combined in various ways, e.g. as mentioned for cases (A)-(E) above.

Acknowledgements

We thank Max Davidson, Simon Grabowsky and Dylan Jayatilaka for extended discussions and for giving us access to the preliminary versions of their manuscripts (Davidson et al., 2022a(Davidson et al., , 2022b) ) before publication.

available from http://journals.iucr.org/services/docxtemplate/) are installed when editing this document.