Hailstorm : A Statically-Typed, Purely Functional Language for IoT Applications

With the growing ubiquity of Internet of Things(IoT), more complex logic is being programmed on resource-constrained IoT devices, almost exclusively using the C programming language. While C provides low-level control over memory, it lacks a number of high-level programming abstractions such as higher-order functions, polymorphism, strong static typing, memory safety, and automatic memory management. We present Hailstorm, a statically-typed, purely functional programming language that attempts to address the above problem. It is a high-level programming language with a strict typing discipline. It supports features like higher-order functions, tail-recursion, and automatic memory management, to program IoT devices in a declarative manner. Applications running on these devices tend to be heavily dominated by I/O. Hailstorm tracks side effects likeI/O in its type system using resource types. This choice allowed us to explore the design of a purely functional standalone language, in an area where it is more common to embed a functional core in an imperative shell. The language borrows the combinators of arrowized FRP, but has discrete-time semantics. The design of the full set of combinators is work in progress, driven by examples. So far, we have evaluated Hailstorm by writing standard examples from the literature (earthquake detection, a railway crossing system and various other clocked systems), and also running examples on the GRiSP embedded systems board, through generation of Erlang.

INTRODUCTION
As the density of IoT devices and diversity in IoT applications continue to increase, both industry and academia are moving towards decentralized system architectures like edge computing [38]. In edge computation, devices such as sensors and client applications are provided greater computational power, rather than pushing the data to a backend cloud service for computation. This results in improved response time and saves network bandwidth and energy consumption [50]. In a growing number of applications such as aeronautics and automated vehicles, the real-time computation is more robust and responsive if the edge devices are compute capable.
In a more traditional centralized architecture, the sensors and actuators have little logic in them; they rather act as data relaying services. In such cases, the firmware on the devices is relatively simple and programmed almost exclusively using the C programming language. However with the growing popularity of edge computation, more complex logic is moving to the edge IoT devices. In such circumstances, programs written using C tend to be verbose, errorprone and unsafe [17,27]. Additionally, IoT applications written in low-level languages are highly prone to security vulnerabilities [7,58].
Hailstorm is a domain-specific language that attempts to address these issues by bringing ideas and abstractions from the functional and reactive programming communities to programming IoT applications. Hailstorm is a pure, statically-typed functional programming language. Unlike impure functional languages like ML and Scheme, Hailstorm restricts arbitrary side-effects and makes dataflow explicit. The purity and static typing features of the language, aside from providing a preliminary static-analysis tool, provide an essential foundation for embedding advanced languagebased security techniques [54] in the future.
The programming model of Hailstorm draws inspiration from the extensive work on Functional Reactive Programming (FRP) [18]. FRP provides an interface to write reactive programs such as graphic animations using (1) continuous time-varying values called Behaviours and (2) discrete values called Events. The original formulation of FRP suffered from a number of shortcomings such as space-leaks [41] and a restrictive form of stream based I/O.
A later FRP formulation, arrowized FRP [46], fixed space leaks, and in more recent work Winograd-Cort et al. introduced a notion of resource types [66] to overcome the shortcomings of the stream based I/O model. The work on resource types is a library in Haskell, and is not suitable to run directly on resource-constrained hardware. Hailstorm uses resource types to uniquely identify each I/O 1 resource. It treats each resource as a signal function to track its lifetime, and prevents resource duplication through the type system. Hailstorm currently has a simple discrete time semantics, though we hope to explore extensions later.
A Hailstorm program is compiled to a dataflow graph, which is executed synchronously. The core of the language is a pure callby-value implementation of the lambda calculus. The synchronous language of arrowized-FRP provides a minimal set of combinators to which the pure core constructs of Hailstorm can be raised. This language of arrows then enforces a purely functional way to interact with I/O, using resource types.
Hailstorm, in its current version, is a work-in-progress compiler which does not address the reliability concerns associated with node and communication failures plaguing edge devices [58]. We discuss the future extensions of the language to tackle both reliability and security concerns in Section 7. We summarize the contributions of Hailstorm as follows: • A statically-typed purely functional language for IoT applications. Hailstorm provides a tailored, purely functional alternative to the current state of programming resource constrained IoT devices. We implement Hailstorm as a standalone compiler, with Erlang and LLVM backends. We have run case studies on the GRiSP embedded system boards [60], to evaluate the features of the language. (Section 4). The compiler implementation and the examples presented in the paper are made publicly available 1 .

LANGUAGE OVERVIEW
In this section we demonstrate the core features and syntax of Hailstorm using running examples. We start with a simple pure function that computes the ℎ Fibonacci number. The simple program above, besides showing the ML-like syntax of Hailstorm, demonstrates some features like (1) higher-order functions (2) recursion (3) partial application (4) tail-call optimization and (5) static typing. All top-level functions in a Hailstorm program have to be annotated with the types of the arguments and 1 https://abhiroop.github.io/ppdp-2020-artifact.zip the return type, which currently allows only monomorphic types. However certain built-in combinators supported by the language are polymorphic which will be discussed in the following section.
The pure core of the language only allows writing pure functions which have no form of interactions with the outside world. To introduce I/O and other side effects, we need to describe the concept of a signal function.

Signal Functions
A fundamental concept underlying the programming model of Hailstorm is that of a Signal Function. Signal Functions, derived from the work on arrowized-FRP [46], are functions that always accepts an input and always returns an output.
Signal functions are analogous to the nodes of a dataflow graph. Signal functions operate on signals which do not have any concrete representation in the language. A signal denotes a discrete value at a give point of time. Nilsson et al [46] use the electric circuit analogy: a signal corresponds to a wire and the current flowing through it, while signal functions correspond to circuit components. An important distinction between Hailstorm and both classic and arrowized-FRP is that signals are always treated as discrete entities in Hailstorm unlike the continuous semantics enforced by FRP.
To create larger programs Hailstorm provides a number of builtin combinators to compose signal functions. These combinators are drawn from the Arrow framework [31] which is a generalization of monads. Arrows allow structuring programs in a point-free style, while providing mathematical laws for composition. We start by presenting some of the core Hailstorm combinators 2 and their types for composing signal functions. Some of the built-in combinators in Hailstorm are polymorphic and the type parameters a, b and c represent the polymorphic types. mapSignal# is the core combinator which lifts a pure Hailstorm function to the synchronous language of arrows, as a signal function (See Fig 1).
Hailstorm then provides the rest of the built-in combinators to compose signal functions while satisfying nine arrow laws [39]. One of the advantages of having a pure functional language is that such laws can be freely used by an optimizing compiler to aggressively inline and produce optimized code. The semantics of composing signals with the arrow combinators is visually depicted in  Now where do signals actually come from? To answer this question -a natural extension to signal functions is using them to interact with I/O, which we discuss in the next section.

I/O
Hailstorm adopts a streaming programming model, where an effectful program is constructed by composing various signal functions in the program. The final program is embedded in a stream of input flowing in and the program transforms that into the output. Response need to be altered every time any new I/O facility has to be added [48]. • There is also no clear mapping between an individual Request and its Response [48].
Work by Winograd-Cort et al. on resource types [66] attempts to address this problem by virtualizing real-world devices as signal functions. What we mean here by "virtualizing" is that the scope of a program is extended so that devices like sensors and actuators are represented using signal functions. For example: sensor : SF () Int uart_rx : SF () Byte actuator : SF Bool () We adopt this model in Hailstorm. The type parameter () represents a void type. There are no values in the language which inhabit the () type. The () type always appears within a signal function. So, an example like sensor : SF () Int -this represents an action which when called, produces an integer.
One of the key aspects of designing a pure functional language is this distinction between an expression that returns a value and an action. When an action (like sensor : SF () Int) is evaluated, it returns a representation of that function call rather than actually executing the call itself. This distinction is the key to equational reasoning in a purely functional program. In the absence of such a distinction the following two expressions are no longer equivalent although they represent the same programs.
--Expression 1 : accepts an input and --duplicates that input to return a pair let x = getInput --makes one I/O call in (x,x) --Expression 2 : accepts two inputs --returns both of them as a pair (getInput, getInput) --makes two I/O calls After enforcing a difference between values and action in the language, we soon encounter one of the pitfalls of treating a realworld object as a virtual device -it allows a programmer to write programs with unclear semantics. For example: def foo : SF () (Int, Int) = sensor &&& sensor Although the above program is currently type correct, it can have two conflicting semantics -(1) either sensor &&& sensor implies two consecutive calls to the sensor device or (2) a single call emitting a pair. Given the type of the sensor function, the latter is not supported and the former is incompatible with Hailstorm's discrete, synchronous semantics.
The notion of Resource Types seeks to solve this problem by labeling each device with a type-level identifier, such that duplicating a device becomes impossible in the program. We change the type of sensor to : The resource keyword in Hailstorm declares a type level identifier which is used for labeling signal functions like sensor above. All the built-in arrow combinators introduced previously are now enriched with new type-level rules for composition as follows: In the combinators above, the type parameters 1 , 2 represent polymorphic resource type variables, which act as labels for the effectful signal functions. The combinator mapSignal# lifts a pure Hailstorm function to a signal function without any effectful operations, and as a result the resource type is Empty. The combinators >>>, &&& and * * * compose signal functions, and result in a disjoint-union of the two resources types. This type-level disjoint union prevents us from copying the same resource using any of these combinators. So in Hailstorm if we try this, def foo : SF ( ∪ ) () (Int, Int) = sensor &&& sensor we currently get the following upon compilation: Type-Checking Error: Error in foo : Cannot compose resources : S S containing same resource Encountered in sensor &&& sensor The type rules associated with composing Hailstorm combinators and their operational semantics are presented formally in Section 3.2 and 3.3 respectively.

Example of performing I/O.
We can distinguish the read and write interface of a resource using two separate resource types. For example, to repeatedly blink an LED we need two APIs -(1) to read its status (2) to write to it. The drivers for these two functions have the following types: The above program runs the function main infinitely. It is possible to adjust the rate at which we want to run this program, discussed later in Section 2.5. This treatment of I/O as signal functions has the limitation that each device (as well as their various APIs) has to be statically encoded as a resource type in the program.

State
Hailstorm supports stateful operations on signals using the loop# combinator. The type of the loop# combinator is slightly different from the type provided by the ArrowLoop typeclass in Haskell, in that it allows initializing the state type variable c. The internal body of the signal function encapsulates a polymorphic state entity. This entity is repeatedly fed back as an additional input, upon completion of a whole step of signal processing by the entire dataflow graph.  The loop# combinator can be used to construct the delay function as found in synchronous languages like Lustre [26], for encoding state.

A sample application
We now demonstrate the use of the Hailstorm combinators in a sample application. The application that we choose is a simplified version of an earthquake detection algorithm [65] which was first used by Mainland et al. to demonstrate their domain specific language for wireless sensor networks [42]. The figure below shows the core dataflow graph of the algorithm. The exponentially weighted moving average (EWMA) component above is a stateful element. We assume the getSample input function is a wrapper around a seismometer providing readings of discrete samples. At the rightmost end, the Detection Event would be another stateful entity which would include some form of an edge detector. The entire program for the earthquake detection is given in Fig 5. The function edge in Fig 5 is a stateful edge detector which generates an action if a boolean signal changes from False to True. To program this, we use an imaginary actuator (like an LED) in a GRiSP board, which would glow red once if the input to it is 1, signalling danger, and otherwise stay green signalling no earthquake.

Sampling rate
The combinators introduced so far execute instantaneously using a logical clock. In Hailstorm, one logical time step includes the following actions, in sequence - • accepting a discrete sample of data from each of its connected input devices • passing the discrete sample through the dataflow graph • finally passing a discrete value to the responsible actuator  A program returning a signal function continuously loops around, streaming in input and executing the above steps, at the speed it takes for the instructions to execute. However under most circumstances we might wish to set a slower rate for the program. The rate# combinator is used for that purpose, The first argument to rate# is the length of the wall clock time (in seconds) at which we wish to set the period of sampling input. This helps us establish a relation between the wall clock time and Hailstorm's logical clock. We demonstrate the utility of the rate# combinator using a Stopwatch example.
2.5.1 Stopwatch. We program a hypothetical stopwatch which accepts an input stream of Ints where 1 represents START, 2 represents RESET and 3 represents STOP. In the above program, the rate# combinator uses the argument g to set the sampling rate to 1 second, which in turn fixes the granularity of the stopwatch as 1 second.

Limitation.
In the current implementation of Hailstorm, an operation like (rate# 2 (rate# 1 sf1)) would result in setting the final sampling rate as 2 , overwriting the value of 1 . In the programs presented here, we use a single clock, and hence a single sampling rate. Libraries like Yampa [14] provide combinators like :: → → which allow oversampling for dealing with multiple discrete sampling rates. As future work, we hope to adopt oversampling operators for communication among signal functions with different sampling rates.

Switches
A Hailstorm dataflow graph allows a form of dynamic, data-driven switching within the graph. It accomplishes this using the switch# combinator: The first argument to switch# accepts a signal function whose output data is used to switch between the various signal functions. The strict typing of Hailstorm restricts the branches of the switch to be of the same type, including the resource type. The switching dataflow is visually presented in  -switch# constrains all of its branches to be of the same type. This is particularly restrictive when dealing with actuators where each actuator would have their own resource type. We currently deal with the notion of choice in the following way -... >>> foo >>> (actuator1 *** actuator2) --instead of emitting a single value foo will emit a pair --of values encoded as (move actuator1, dont move actuator2) The Arrow framework provides more useful combinators based on the ArrowChoice typeclass which are currently absent from Hailstorm.
-The switch# combinator is an experiment to describe expressions of the form: The combinator currently supports expression only of the above form. Thus, we currently do not allow more general functions of type (b→SF b c) as the second parameter to switch#, and so avoid problems with possibly undefined runtime behaviour. However, there is no check in the compiler to enforce this restriction. As future work we hope to adopt the guards syntax of Haskell to represent the second parameter as a collection of boolean clauses and their corresponding actions.

SYNTAX, SEMANTICS AND TYPES
In the previous sections, we have given an informal treatment to the most important syntactic parts of Hailstorm, for IoT applications, and described the programming model using them. In this section, we present the core syntax, type rules and operational semantics of the main parts of Hailstorm.

Syntax
The set of types in the source language is given by the following grammar. : The type 1 =⇒ 2 represents a signal function type from a to b with the resource type r i.e SF r a b.
The abstract syntax of the core expressions of Hailstorm is given by the following grammar. The meta-variable ∈ ranges over variables of the source language. Additionally we let ∈ Z and ∈ R. We use and separately to denote ordinary expressions and arrow based signal function expressions respectively.
In the grammar above we describe two primitives for I/O called read# and write#. In practise, as Hailstorm deals with a number of I/O drivers there exists a variety of I/O primitives with varied parameters and return types. However for the purpose of presenting the operational semantics, we abstract away the complexity of the drivers and use the abstracted read# and write# to describe the semantics in Section 3.3.

Type rules
Hailstorm uses a fairly standard set of type rules except for the notion of resource types. The typing context of Hailstorm employs dual contexts, in that it maintains (1) Γ -a finite map from variables to their types and (2) Δ -a finite set which tracks all the I/O resources connected to the program. Δ; Γ ::= · | Γ, : An empty context is given by ·. Additionally (Γ) provides the set of variables bound by a typing context.
In Fig 7, we show the most relevant type rules concerning signal functions and their composition. The remaining expressions follow standard set of type rules which is provided in its entirety in the extended version of this paper [55].
Looking at the rule T-Mapsignal, a signal function such as 1 ∅ =⇒ 2 denotes the result of applying mapSignal# to a pure function. This results in an expression with an empty resource type denoted by ∅. A missing rule is the introduction of a new resource type in the resource type context Δ. The resource type context is an appendonly store and a new resource is introduced using the keyword resource. It can be defined using this simple reduction semantics where denotes one step of reduction which occurs at compiletime. The rules such as T-Compose, T-Fanout, T-Combine, T-Switch apply a type-level disjoint union to prevent resource duplication.

Big-step Operational Semantics
In this section, we provide a big-step operational semantics of our implementation of the Hailstorm language, by mapping the meaning of the terms to the lambda calculus. We begin by defining the values in lambda calculus that cannot be further reduced: i and f represent integer and float constants respectively. We use n to represent variables and the last term denotes a lambda expression. The syntax that we use for defining our judgements is of the form : The variables 1 , 2 are finite partial functions from variables n to their bound values ∈ . In case a variable n is unbound and s is called with that argument it returns ∅. The above judgement is read as starting at state 1 and evaluating the term M results in the irreducible value ∈ while setting the final state to 2 .
The first judgement essential for our semantics is this, which means that the values V cannot be reduced further. We use a shorthand notation ( , ) to signify lookup the variable n in s. Additionally, to make the semantics more compact we  They do not belong to V but it is possible to represent all three of them using plain lambdas and function application -shown in the extended version of this paper [55].
In Fig 8, we show the most relevant big-step operational semantics concerning signal function based combinators. The remaining expressions have standard semantics and the complete rule set is provided in the extended paper [55]. In the rule eval-Rate, we use Ψ to store the sampling rate. In our current implementation, when composing signal functions with different sampling rates, the state transition from 2 to 3 overwrites the first sampling rate.
In eval-Loop-Init and eval-Loop, the subscript n represents a variable name that is used as a key, in the global state map s, to identify each individual state. Each Hailstorm program compulsorily has a main function. After the entire program is partially evaluated (described in Section 4.1) and given that the main function causes a side effect (denoted by () below), we have the following rule: The eval-Main rule demonstrates the streaming and infinite nature of a Hailstorm program when the main function is a signal function itself. After causing a side effect, it calls itself again and continues the stream of effects while evaluating the program using the semantics of Fig 8.

IMPLEMENTATION
Here we describe an implementation of the Hailstorm language and programming model presented in the previous sections. We implement the language as a compiler -the Hailstorm compilerwhose compilation architecture is described below. The compiler pipeline starts by parsing a Hailstorm source file and desugaring small syntactic conveniences provided to reduce code size (such as case expressions). After desugaring, the Hailstorm AST constituted of the grammar described in Section 3.1 is generated. Next, typechecking of the AST using the type rules of Fig 7 is done. After a program is typechecked, it is transformed into the Hailstorm core language called L1.
The L1 language is an enriched version of the simply typed lambda calculus (STLC) with only nine constructors. Unlike STLC, L1 supports recursion by calling the name of a global function. It is currently incapable of recursion using a let binding. L1 has a simpler type system than the Hailstorm source, as it erases the notion of resource types -which are exclusively used during type checking.
The L1 language being considerably simpler, forms a sufficient foundation for running correctness preserving optimization passes. L1 attempts to get rid of all closures with free variables, as they are primarily responsible for dynamic memory allocation. Additionally, it attempts to inline expressions (primarily partial applications) to further reduce both stack and heap memory allocation. The optimization passes are described in further detail in Section 4.2.
The final optimization pass in L1 is a partial evaluation pass which specializes the program to convert it into an effectful program. Before this pass, all I/O inducing functions are treated as values, by guarding them inside a abstraction. This pass evaluates those expressions by passing a compile-time token and turning the L1 non-effectful program into an effectful one. This pass is discussed in further detail in Section 4.1.
Finally, the effectful Hailstorm program gets connected to one of its backends. Recursion is individually handled in each of the backends by using a global symbol table. We currently support an LLVM backend [36] and a BEAM backend [4].

I/O
The I/O handling mechanism of the Hailstorm compiler is the essential component in making the language pure. A pure functional language allows a programmer to equationally reason about their code. The entire program is written as an order independent set of equations. However, to perform I/O in a programming language, it is necessary to (1) enforce an order on the I/O interactions, as they involve chronological effects visible to the user and (2) interact with the real world and actually perform an effect.
To solve (1) we use the eval-App rule given in Fig 8. Hailstorm, being a call-by-value implementation of the lambda calculus, follows ordering in function application by always evaluating the function argument before passing it to the function. This semantics of function application allows us to introduce some form of ordering to the equations.
To solve (2) we have extended our pure lambda calculus core to involve effectful operations like read and write. Again, as observed in the semantic rules, eval-Read and eval-Write from Fig 8, the effectful operations are guarded by abstractions to treat them as values, rather than operations causing side-effects. This allows us to freely inline or apply any other optimization passes while preserving the correctness of the code. However, to finally perform the side effect, we need to resolve this lambda abstraction at compile time. We do this by partially evaluating [33] our program. Below we describe the semantics of the partial evaluation step.
A Hailstorm program is always enforced (by the typechecker in our implementation) to contain a main function. We shall address the case of an effectful program i.e a main function whose return type is an effectful signal function such as SF r () (). The () type in the input and output parameters reflect that this function reads from a real world source like a sensor and causes an effect such as moving an actuator. The partial application pass is only fired when the program contains the () type in either one of its signal function parameter.
Let us name the main function above with the return type of SF r () () as . Now is itself a function, embedded in a stream of input and transforming the input to an output, causing an effect. As it is a function we can write: = . . Additionally in our implementation, after typechecking, a signal function type 1 =⇒ 2 is reduced to a plain arrow type -1 → 2 . Now, using the two aforementioned definitions, we can write the following reduction semantics: = . : where denotes one step of reduction and denotes an arbitrary compile time token of type (). The final step, which produces the expression [ ↦ → : ()], is the partial evaluation step. We demonstrate this reduction semantics in action using an example: Given the above definition of main, we can apply the reduction rule eval-Partial and eval-App to get the following, (Partial evaluation) Now the program is ready to create a side effect as the read function is no longer guarded by a abstraction. The eval-App rule guarantees that read is evaluated first and only then the value is fed to write owing to the call-by-value semantics of Hailstorm.
4.1.1 Limitation. The Hailstorm type system doesn't prevent a programmer from writing write# >>> read#. The type of such a program would be SF (STDOUT U STDIN) Int Int. As the types do not reflect the () type, the partial evaluation pass is not fired and the program simply generates an unevaluated closurewhich is an expected result given the meaningless nature of the program. However, the type would allow composing it with other pure functions and producing bad behaviour if those programs are connected to meaningful I/O functions. This is simply solved by adding the following type rule: Our implementation treats HOFs as closures which are capable of capturing an environment by allocating the environment on the heap. In line no. 3 above, the value of the variable x is heap allocated. However, in resource constrained devices heap memory allocation, is highly restrictive and any language targeting such devices should attempt to minimize allocation. Hailstorm applies a lambda-lifting [32] transformation to address this. Lambda-lifting lifts a lambda expression with free variables to a top-level function and then updates related call sites with a call to the top-level function. The free variables then act as arguments to the function. This effectively allocates them on the stack (or registers). Owing to our restricted language and the lack of polymorphism, our algorithm is less sophisticated than the original algorithm devised by Johnsson. We describe the operational semantics of our algorithm in the lambda-lifting rule below.
We use a slightly different notation from Section 3.3 here. is used to describe the entire program with its collection of top level functions. We identify a modified program with [ ( ). ] to mean a new program with an additional global function G which accepts an argument x and returns an expression E. Finally the function is used find the set of free variables in a Hailstorm expression. ..., , ) . ] The above rule returns a modified expression ′ which consists of a function call to ( 1 , ..., , ) . with the free variables 1 , ..., as arguments. This rule is repeatedly run on local lambda expressions until none has any free variables. It transforms the program foo above to : Inlining subsumes optimization passes like copy-propagation in a functional language. Our prototype inliner is relatively conservative in that (1) it doesn't attempt to inline recursive functions and (2) it doesn't attempt inter-function inlining.
However, it attempts to cooperate with the lambda-lifter to remove all possible sites of partial applications, which are also heap allocated, to minimize memory allocation. The program shown in the previous section, undergoes the cycle in Fig 10   wrapper around the driver for a simulated seismometer, which when polled for data provides a reading. The rightmost edge of the graph for the Detection Event pulls on the dataflow graph after it completes an action and the rest of the graph in turn pulls on getSample, which polls the simulated seismometer.
However, in certain devices such as UART, a push based model is more prevalent, where data is asynchronously pushed to the drivers. In such cases, to avoid dropping data the wrapper function (such as getSample) needs to be stateful and introduce buffers that store the data. The pseudocode of the wrapper function for such a driver, written for our Erlang backend, would look like Fig 12. loop (

Figure 12: Enforcing pull semantics on push-based data
In Fig 12 there are two separate message calls handled. The data transmission from the drivers is handled using the uartdriver message call, which continuously buffers the data. On the other hand, the Hailstorm program, upon finishing one cycle of computation, requests more data using the hailstormcall message, and proceeds with the rest of the cycle.

Limitation.
Elliott has criticized the use of pull semantics [19] as being wasteful in terms of the re-computation required in the dataflow graph. He advocates a hybrid push-pull approach, which, in case of continuously changing data, adopts the pull model, but in the absence of any change in the data doesn't trigger the dataflow graph. This approach could be useful in resource constrained devices, where energy consumption is an important parameter, and we hope to experiment with this approach in future work.

The digital -analog interface
The runtime of Hailstorm has to deal with the boundary of discrete digital systems and continuous analog devices. The input drivers have to frequently discretize events that occur at some unknown point of time into a stream of discrete data. An example is the Stopwatch simulation from Section 2.5.1. The pressing of an ON button in a stopwatch translates to the stream of ones (111...) and when switched OFF the simulation treats that as a stream of zeroes. This stream transformation is handled by the wrapper functions around the input drivers.
On the contrary, for the output drivers a reverse translation of discrete to continuous is necessary. We can take the example of operating traffic lights (related demonstration in Section 5.1.2). When operating any particular signal like GREEN supplying a discrete stream of data (even at the lowest granularity) will lead to a flickering quality of the light. In that case the wrapper function for the light drivers employs a stateful edge detector, as discussed in Section 2.4, to supply a new signal only in case of change.

Backend specific implementation
The backend implementation includes • Memory management. The compiler attempts to minimize the amount of dynamically allocated memory using lambdalifting and inlining such that the respective garbage collectors have to work less. Future work hopes to experiment with static memory management schemes like regions [61]. • Tail call optimization (TCO). Erlang itself does TCO and LLVM supports TCO when using the fastcc calling convention.

EVALUATION 5.1 Case Studies
In this section we demonstrate examples from the synchronous language literature [30] written in Hailstorm.

Watchdog process.
A watchdog process monitors a sequential order processing system. It raises an alarm if processing an order takes more than a threshold time. It has two input signals -(1) : () which emits when an order is placed and otherwise, (2) : () which also emits only when an order is done. For output we use -: () where an alarm is rung only when is supplied. In the program below we keep a threshold time for order processing as 3 seconds. 5.1.2 A simplified traffic light system. We take the classic example of a simplified traffic light system from the Lustre literature [53].
The system consists of two traffic lights, governing a junction of two (one-way) streets. In the default case, traffic light 1 is green, traffic light 2 is red. When a car is detected at traffic light 2, the system switches traffic light 1 to red, light 2 to green, waits for 20 seconds, and then switches back to the default situation. We use a sensor -: () -which keeps returning as long as it detects a car. The system, upon detecting a subsequent car, resets the wait time to another 20 seconds. We sample from the sensor every second. For the traffic lights, we use 1 to indicate green and 0 for red. The problem. We consider a two-track railway level crossing area that is protected by barriers, that must be closed in time on the arrival of a train, on either track. They remain closed until all trains have left the area. The barriers must be closed 30 seconds before the expected time of arrival of a train. When the area becomes free, barriers could be opened, but it's not secure to open them for less than 15 s. So the controller must be warned 45 s before a train arrives. Since the speed of trains may be very different, this speed has to be measured, by detecting the train at two points separated by a known distance. A first detector is placed 2500 m before the crossing, and a second one 100 m after this first. A third is placed after the crossing area, and records a train's leaving. We divide our solution into three programs.
The Detect Process. The passage of a train is detected by a mechanical device, producing a pulse -: () only when a wheel runs on it (otherwise ). The detect process receives all these pulses, but warns the controller only once, on the first wheel. All following pulses are ignored.  Given the alarms from the two separate tracks, the barrier controller sends an open/close signal represented by 0 and 1 respectively. In case a train is approaching in both of the tracks at the same speed -the barrier for only track 1 is opened. Note that the three instances of rate# all sample at the same interval of 0.1 seconds. Hailstorm currently doesn't have well defined semantics for programs with multiple clock rates.

Microbenchmarks
Here we provide memory consumption and response-time microbenchmarks for the examples presented above using the Erlang backend.   Table 1 shows the response time for the programs measured in milliseconds. We measure the CPU Kernel time (CPUT) -which calculates the time taken by the dataflow graph to finish one cycle of computation. We show three separate runs where the Erlang virtual machine is killed and restarted to reset the garbage collector. Each of the numbers are an average of twenty iterations of data processing. We use the erlang:statistics module for measuring time and in the applications I/O happens over the command line interface, which explains the overall slow behaviour (tens of milliseconds). The metrics are run on the erts-10.6.4 runtime and virtual machine running on a Macbook-Pro with a 2.9 GHz Intel Core i9 processor. A common observation is that the computation takes less than 1% of the total wall clock time involved in the response rate, showing that the I/O reading/writing times dominate the final response rate.
An alternate benchmarking strategy which we used was to model the input from the sensors as an in-memory structure (in Erlang) and compute the total response time for processing those values using the timer:tc module: The above values are all measured in microseconds which are averaged over forty iterations each. As expected, using an in-memory structure results in graph processing times that are much lower than those in Table 1.

RELATED WORK 6.1 Programming Languages for IoT
There has been recent work on designing embedded DSLs (EDSLs) for IoT applications [11]. In EDSLs, I/O is handled by embedding a pure core language inside a host language's I/O model -which is an approach that Hailstorm deliberately avoids. Given the I/O dominated nature of IoT apps, we choose to focus much of our attention on designing a composable stream based I/O model, rather than only considering the pure core language, as many EDSLs do.
Other approaches like Velox VM [62] runs general purpose languages like Scheme on specialized virtual machines for IoT devices. A separate line of work has been exploring restrictive, Turingincomplete, rule-based languages like IoTDSL [2] and CyprIOT [9].
Juniper [29] is one of the few dedicated languages for IoT but it exclusively targets Arduino boards. Emfrp [56] and its successor XFRP [59] are most closely related to the goals of Hailstorm. However, their model of I/O involves writing glue code in C/C++ and embedding the pure functional language inside it. Hailstorm has a more sophisticated I/O integration in the language.
While IoT stands for an umbrella term for a large collection of software areas, there has been research on declarative languages for older and specialized application areas like: • Wireless Sensor Networks (WSNs). There exists EDSLs like Flask [42] and macroprogramming languages like Regiment [45] and Kairos [24] for WSNs. • Real Time Systems. Synchronous language like Esterel [10], Lustre [26] are a restrictive set of languages designed specifically for real time systems. Further extensions of these languages like Lucid Synchrone [12], ReactiveML [44], Lucy-n [43] have attempted to makes them more expressive. The applications demonstrated in the paper are expressible in synchronous languages, albeit using a very different interface from Hailstorm. While languages like Lustre and its extensions are pure they restrict their synchronous calculus to the pure core language and handle I/O using the old stream based I/O model of Haskell [48]. Hailstorm explores the design space of pure functional programming with the programming model and purity encompassing the I/O parts as well.
In Lustre, a type system called the clock calculus ensures that programs can run without any implicit buffering inside the program. Strong safety properties such as determinism and absence of deadlock are ensured at compile time, and programs are compiled into statically scheduled executable code. This comes at the price of reduced flexibility compared to synchronous dataflow-like systems, particularly in the ease with which bounded buffers can be introduced and used. Mandel et al have studied n-synchronous systems in an attempt to bring greater flexibility to synchronous languages [43]. Hailstorm, in the presence of recursion, is unable to statically predict memory usage but we plan future work on type level encoding of buffer sizes to make memory usage more predictable. Polychronous languages like SIGNAL [8] and FRP libraries like Rhine [6] provide static guarantees on correctness of systems with multiple clocks -something that we hope to experiment with in the future.
Various implementations have explored the choice between a static structure of the dataflow graph (for example Elm [15]) or dynamic structure, as in most Haskell FRP libraries [3] as well as FrTime in Racket [13]. The dynamic graph structure makes the language/library more expressive, allowing programs like sieves [25].
The higher-order FRP implementation offers almost the local maxima of tradeoffs, but at the cost of an extremely sophisticated type system, which infests into the source language, compromising its simplicity.
The loop# combinator in Hailstorm is similar to the -combinator first introduced by Sheeran [57] and to the loopB combinator in the Causal Commutative Arrows (CCA) paper [40]. CCA presents a number of mathematical laws on arrows and utilizes them to compile away intermediate structures and generate efficient FRP code -a promising avenue for future work in Hailstorm.
FRP has seen adoption in various application areas. Application areas related to our target areas of IoT applications include robotics [64], real-time systems [63], vehicle simulations [21] and DSLs discussed in detail in the previous section.

Functional I/O
A detailed summary of approaches to functional I/O was presented by Gordon et al. [22]. Since then monadic I/O [23] has become the standard norm for I/O in pure functional languages, with the exception of Clean's I/O system [1] based on uniqueness types.
More recently, there have been attempts at non mondaic I/O using a state passing trick in the Universe framework [20]. The latest work has been on the notion of resource types, proposed by Winograd-Cort et al. [66], which is explored as a library in Haskell. In effect, their library uses Haskell's monadic I/O model; however, they mention possible future work on designing a dedicated language for resource types. Hailstorm explores that possibility by integrating the idea of resource types natively in a language's I/O model. FRPNow! [49] provides an alternate monadic approach to integrate I/O with FRP.

FUTURE WORK
Hailstorm is an ongoing work to run a pure, statically-typed, functional programming language on memory constrained edge devices. As such, there are a number of open avenues for research: -Security. We hope to integrate support for Information Flow Control (IFC) [28] -which uses language based privacy policies to determine safe dataflow -in Hailstorm. Given that the effectful combinators are based on the Arrow framework, we expect to build on two lines of work -(1) integrating IFC with the Arrow typeclass in Haskell [37] (2) a typeclass based technique to integrate IFC without modifying the compiler [52] but one which relies on the purity and static typing of the language.
-Reliability. Hailstorm currently doesn't provide any fault tolerance strategies to mitigate various node/communication failures. However, being hosted on the Erlang backend, we plan to experiment with distributed versions of the arrow combinators where the underlying runtime would use supervision trees to handle failures. This would be a much more invasive change as the synchronous dataflow model of the language is not practically suitable in a distributed scenario -leading to more interesting research directions on macro-programming models [24].
-Memory constrained devices. The evaluation of this paper is carried out on GRiSP boards which are relatively powerful boards. We are currently working on developing a small virtual machine which can interpret a functional bytecode instruction set and run on much more memory constrained devices like STM32 microcontrollers.

CONCLUSION
We have presented the design and implementation of Hailstorm, a domain specific language targeting IoT applications. Our evaluation suggests that Hailstorm can be used to declaratively program moderately complex applications in a concise and safe manner. In the future, we hope to use the purity and type system of Hailstorm to enforce language-based security constraints, as well as to increase its expressiveness to enable the description of interacting IoT devices in a distributed system.


extfieldsofstudy:
Computer Science
Physics
Mathematics
Medicine

text:
Molecule Edit Graph Attention Network: Modeling Chemical Reactions as Sequences of Graph Edits

The central challenge in automated synthesis planning is to be able to generate and predict outcomes of a diverse set of chemical reactions. In particular, in many cases, the most likely synthesis pathway cannot be applied due to additional constraints, which requires proposing alternative chemical reactions. With this in mind, we present Molecule Edit Graph Attention Network (MEGAN), an end-to-end encoder-decoder neural model. MEGAN is inspired by models that express a chemical reaction as a sequence of graph edits, akin to the arrow pushing formalism. We extend this model to retrosynthesis prediction (predicting substrates given the product of a chemical reaction) and scale it up to large datasets. We argue that representing the reaction as a sequence of edits enables MEGAN to efficiently explore the space of plausible chemical reactions, maintaining the flexibility of modeling the reaction in an end-to-end fashion, and achieving state-of-the-art accuracy in standard benchmarks. Code and trained models are made available online at https://github.com/molecule-one/megan.

Introduction
Synthesis planning answers the question of how to make a given molecule. Due to the substantial size and complexity of the reaction space, synthesis planning is a demanding task even for skilled chemists and remains an important roadblock in the drug discovery process 1 . Computer-aided synthesis planning (CASP) methods can speed-up the drug discovery process by assisting chemists in designing syntheses [2][3][4][5] . These methods are also increasingly employed in de novo drug design, in which estimating the synthetic accessibility of a large number of compounds has surfaced as a key challenge to the field 6 .
Designing a synthesis plan involves predicting reaction outcomes for a given set of possible substrates (forward synthesis prediction), or proposing reactions that can simplify a given target molecule (retrosynthesis prediction) 7 . In both cases, it is key for automated synthesis planning tools to model reliably outcomes of a diverse set of chemical reactions. This is especially important when there are additional constraints imposed on the desired synthesis plan such as avoiding certain starting materials or using green chemistry 8 .
Many approaches to predicting reaction outcomes employ a static library of reaction templates.
A reaction template encodes a graph transformation rule that can be used to generate the reaction based on the product or subtrates [9][10][11][12][13][14] . Due to computational limitations, they typically require representing reactions using a restricted number of templates, which necessarily limits the coverage of the chemical reaction space accessible by such methods.
The small coverage and other limitations of template-based methods motivated the development of more flexible deep learning models, which aim to learn transformations and their applicability rules directly from the data 15,16 . Despite the added flexibility and the resulting better coverage of the chemical space, they tend to propose a limited set of plausible chemical reactions for a given input 17 .
Models such as proposed in Schwaller et al. 15 , Karpov et al. 16 belong to the class of sequence to sequence models. They generate a chemical reaction by sequentially outputting individual symbols in the SMILES notation 18 . We hypothesize this makes it challenging to generate more diverse predictions from the model, as there is no natural decomposition of the predictive distribution into plausible but different reactions. Notably, many previous papers have pointed out limited diversity of sequence to sequence models in other domains [19][20][21] .
An arguably more natural idea is to express a reaction as a sequence of graph edits such as bond addition or removal, which is inspired by how chemists describe reactions using the arrow pushing mechanism 22,23 . Bradshaw et al. 22 introduce a model that predicts the reaction mechanism represented as a sequence of bond removal and additions that directly correspond to electron paths.
However, their approach is limited to a certain subset of the chemical reaction space (reactions with linear chain topology) and forward synthesis. In a similar spirit, Do et al. 23 models a reaction as a set of operation on atom pairs. However, due to the lack of support for atom addition, their method also cannot be readily applied to retrosynthesis prediction. We hypothesize that expressing reaction in this way enables to easily search through different alternative reactions because the proposed reaction is obvious usually after the first few generation steps. This allows the generation algorithm to naturally backtrack and explore other possible reactions. We note that related methods have a rich history in parallel applications of generative deep learning such as natural language processing 24 .
In this work, we propose the Molecule Edit Graph Attention Network (MEGAN), an encoderdecoder end-to-end model that generates a reaction as a sequence of graph edits. Designing a novel architecture and training procedure enables us to apply this approach to both retro and forward synthesis, and scale it up to large datasets. More specifically, our two main contributions are as follows: • We extend the idea to model the reaction as a sequence edits to retrosynthesis prediction and scale it to large datasets. This required introducing a novel encoder-decoder architecture and an efficient training procedure that avoids the need to use reinforcement learning.
• We achieve competitive performance on retrosynthesis and competitive performance on forward synthesis as well as state-of-the-art top-k accuracy for large K values on all tested datasets. This serves as evidence that MEGAN achieves excellent coverage of the reaction space.
We open-sourced our code and trained models at https://github.com/molecule-one/megan.

Related work
Computer-aided reaction prediction has a rich history 25 . Early approaches to generating reactions relied on manually crafted rules 26 , which were difficult to apply to novel chemistry. Another family of methods uses physical chemistry calculations 27 . These methods are typically too computationally intensive to be used in synthesis planning software and can be complemented using statistical methods that learn from data.
Statistical approaches to predicting reactions and designing synthesis paths can be broadly categorized into template-based and template-free approaches. Template-based approaches use reaction rules or templates. These templates can be automatically mined from a database of known reactions or defined manually [9][10][11][12][13][14] .
Manually encoding reaction templates necessitates developing complex rules governing their applicability 14,28 . This motivated the recent rise in popularity of machine-learning models used to select or rank most relevant templates 3,29 . Such models can be applied to assess the reactivity of atoms to which reaction rules should be applied 30,31 . An interesting alternative to this approach is to prioritize transformations that have already been applied to similar molecules 32 . A suitable set of templates combined with a strong ranking model can achieve state-of-the-art performance on standard retrosynthesis benchmarks. 13 The small coverage and other limitations of template-based methods motivated the development of template-free models. A particularly successful class of template-free methods is based on sequence to sequence models that sequentially (symbol by symbol) predict the target SMILES 18 string from the input SMILES string 18 . Such approaches can be used out-of-the-box for both forward synthesis 15,33 and retrosynthesis 16,34 . However, they have been shown to make some trivial mistakes 28 , and produce reactions with a limited diversity 17 . These models also act as black-boxes; they do not provide the reasoning behind their predictions and are not able to map atoms between the substrates and the product.
We hypothesize that the limited diversity of sequence to sequence models can be attributed to the fact that they generate reaction as a sequence of SMILES symbols from left to right. Many papers have also pointed out shortcomings in the diversity of the generated output by sequence to sequence models [19][20][21] . In natural language processing, there is a rich parallel history of exploring alternative generation ordering and adding through various means structure to the output. Gu et al. 24 show that using an alternative order of output symbols than left to right enables generating more diverse prediction. Mehri and Sigal 35 propose an alternative decoding algorithm to improve the diversity of sequence to sequence models that starts from generating the middle word. These modeling choices can be viewed as incorporating more suitable inductive biases for the given domain. From this perspective, we study an alternative representation of the output in generative models applied to reaction outcome prediction.
Our work focuses on a recently introduced class of template-free methods, which defines reaction generation as predicting target graph by sequentially modifying the input graph. Such an approach, while remaining template-free, can provide greater interpretability of predictions, which are modeled as direct transformations on molecules. Models that predict graph edits have been successfully applied to forward synthesis 23 and can be fine-tuned to provide especially interpretable results on certain subsets of the chemical reaction space 22 . In this work, we propose a model inspired by these prior works and extend the approach to retrosynthesis prediction.
Concurrently to our work, the graph-edit reaction generation approach has been employed to retrosynthesis in G2GS and GRAPHRETRO models [36][37][38] . The G2GS model presented in Shi et al. 36 , Yan et al. 38 consists of two separate modules for predicting the reaction center and generating final substrates from the disconnected synthons. GRAPHRETRO 37 employs a similar framework but completes the synthons with substructures which are selected from a predefined set found on the training data. In contrast to these concurrent works, we propose an end-to-end model for generating reaction through a sequence of edits and apply it both to retro and forward synthesis prediction.
Our work on reaction outcome prediction is motivated by their wide utility. An important use of reaction outcome prediction models arises in automated synthesis planning; see for example Segler et al. 3 for an excellent discussion. It is also worth noting that reaction outcome prediction has recently risen in importance in the context of de novo drug discovery. Gao and Coley 39 have drawn attention to the fact that many recent de novo drug discovery methods tend to generate hard to synthesize molecules. In this context, automated synthesis planning can be used to filter out these hard to synthesize compounds, and compares favorable to other simpler methods 40,41 . Other works use reaction outcome prediction models to define the space of accessible molecules to de novo models 42 .

Molecule Edit Graph Attention Network
Molecule Edit Graph Attention Network (MEGAN) is an encoder-decoder architecture based on the graph convolutional network architecture that predicts a sequence of actions on atoms and bonds of a chemical compound.

Input and output representation
We begin by describing the input and the output representation of the MEGAN model.  We add an additional node to the graph named supernode 44 , which is connected to all atoms in the graph with a special SUPERNODE bond type. In graph convolutional neural networks, supernode is particularly useful for passing information between connected components of a graph, especially after it was split by deleting a bond. Adding the supernode makes n equal to the number of atoms in the input molecule + 1.
To construct feature vectors for atoms and bonds , we use the following chemical properties calculated using the RdKit package 45 : ATOMIC NUMBER, FORMAL CHARGE, CHIRAL TAG, NUMBER OF EXPLICIT HYDROGEN ATOMS, IS AROMATIC for atoms and BOND TYPE, BOND STEREOMETRY (None, Z or E) for bonds. All these features are set to zero for the supernode and all bonds with the supernode, except for BOND TYPE which has a special supernode value for bonds between atoms and the supernode. We also add a binary flag to indicate whether a given atom or bond was edited or added in any of the previous steps. Finally, we indicate whether a given node is an atom or the supernode with a binary flag on each node. Final atom and bond vectors are constructed by concatenating one-hot vectors of all of their categorical features. We describe the remaining details in Section Featurization in the Supplement.
Output representation The output of MEGAN is a sequence of graph actions, which are then applied to the input graph to produce the desired output graph. An action can be applied to a single atom or a pair of atoms. In the latter case the action modifies the bond between the two atoms.
To train the model, we compare the outputted sequence of actions with the ground-truth sequence based on a mapping between the product and the substrates (see Section Gradient-based training of MEGAN). We define the following graph actions: • Edit atom properties (EditAtom) • Edit bond between two atoms (EditBond) • Add new atom to the graph (AddAtom) • Add new benzene ring to the graph (AddBenzene)

• Stop generation (Stop)
EditAtom changes properties of atoms, such as the formal charge, chirality, or aromaticity.
EditBond adds, edits, or deletes a bond between two atoms. AddAtom adds a new atom of a specified type as a neighbor of another atom already existing in the graph, with a specified bond type. AddBenzene simplifies generation by allowing the model to append a complete benzene ring to a selected carbon atom. Stop action indicates the end of the generation process and can be outputted only by the supernode. We describe the possible actions in detail in the "Graph edit actions" in Supplement. f

Model architecture
Encoder-decoder architecture of MEGAN The MEGAN model consists of two parts: the encoder, which is invoked only once per reaction generation at its beginning, and the decoder, which is sequentially invoked to generate actions. This general encoder-decoder architecture is shared by the most successful sequence to sequence models in deep learning Vaswani et al. 46 .
Both encoder and decoder are variants of graph convolutional networks (see GCN-att layer) that modify hidden vector representation of nodes. They take the same input and output node feature vector size n a and input edge feature vector size n b . The input graph x t is modified according to the action predicted by the decoder at each step t, until a special Stop action is outputted. We stack n e and n d layers of the encoder and decoder, respectively. The architecture and the inference process are depicted schematically in Figure 1.
Generating graph actions At each generation step t, the one-hot encoded features of the current A t = g emb (A OH t ) to convert atoms and bonds to n a and n b -dimensional spaces, respectively. We treat n a and n b as hyperparameters.
In the first generation step (t = 1), we pass the embedded graph through the encoder and the decoder to arrive at atom features H 1 = dec(enc(H 0 1 )). At each subsequent step t, we only use the decoder. We incorporate previous step representation by taking the element-wise maximum operation over the current and previous atom features. Namely, for t > 1, H t = dec(max(enc(H t−1 ), H t−1 )).
We zero-pad features of H t−1 for any node that was added to the graph at step t. The encoder and decoder also take as input A t .
The logits L a t and L b t for atom actions Act a and bond actions Act b are calculated at each step t from H t and A t as follows: , where g atom , g atom , g bond , and g bond are linear layers; σ r is the ReLU activation function. The GCN-att layer The basic building block of the encoder and decoder is an attention-based GCN layer that we call GCN-att. We enhance the GCN layer from Veličković et al. 43 by adding bond features as input information for computing the attention values. Let H t ∈ R n×h denote input node features for the t-th GCN-att layer and N (i) ⊂ Z ≥0 denote set of indices of neighbors of node at index i (where i ∈ N (i)). We calculate new node features H t+1 ∈ R n×h as follows: where as before σ r denotes ReLU activation function and indicates vector concatenation. The scalar K ∈ N + is the number of attention heads and f att : f : R na → R na/K are standard linear layers. Numbers n a , n b , d and K are hyperparameters of the model. We require that n a is divisible by K. We use the same hyperparameter values for all GCN-att layers in the model.

Gradient-based training of MEGAN
In contrast to Do et al. 23 , who use reinforcement learning to train their model, we back-propagate directly through the maximum likelihood objective to train MEGAN. We use teacher forcing to train the model, i.e. to predict each step during reaction generation we use previous steps from the ground-truth as input to the model. This is nontrivial, as computing the gradient of the likelihood objective requires defining a fixed ordering of actions 47 . To solve this issue, You et al. 47 enumerates atoms using breadth-first search.
We adapt a similar idea to reaction generation.
At first, we define a general ordering of action types for forward-and retrosynthesis. This ordering defines which types of actions should be performed first if there are more than one type of actions that could be applied to an atom in a generation step. There are various possible orderings, described in detail in the Section "Ablation study on action ordering". The highest accuracy is achieved in the ordering called BFS RAND-AT. In this ordering, for retrosynthesis, we give bond deletion the highest priority, as it is usually the step that determines the reaction center, which is the first step of reaction prediction in other methods. [36][37][38]48 Analogously, adding a bond has the highest priority for forward synthesis, as this usually determines the reaction center in forward prediction.
Priorities of other types of actions were determined experimentally; we did not observe a significant change in validation error when modifying these priorities. We list these priorities explicitly in Table 1.
The remaining ties between actions are broken as follows. For retrosynthesis prediction, we prioritize actions that act on the least recently modified atoms. Our motivation is to prioritize more difficult actions (such as deciding where the reaction occurs) at the expense of simpler actions (such as completing the synthons resulting from breaking a bond). In forward prediction, we prioritize the most frequently visited atom, which we found to perform better empirically. The remaining ties are broken at random. We describe the algorithm in full detail in the Supplement (Ground truth action ordering).
It is worth noting that the action ordering algorithm is used only to acquire ground truth samples for training MEGAN with a gradient-based method. The model is expected to learn a strong prior towards outputting edits in such order; however, when evaluating MEGAN, we do not explicitly limit the available action space, giving it a chance to determine the order by itself. We observe that MEGAN often predicts the ground truth reaction using a different order of actions than the ground truth training sample.
We use the following rules to resolve conflicts when the action could be outputted at different atoms or bonds: • Stop action can be predicted only by the supernode • Bond actions can be predicted for indices i and j, where i < j and nodes at i and j are atoms Finally, for simplicity, we do not mask out redundant actions, such as deleting a non-existing bond or editing the atom to the same values of properties, as we expect the model to learn not to use such actions.

Experiments
We evaluate MEGAN on three standard datasets for retro and forward synthesis prediction. We first evaluate it on the standard retrosynthesis prediction benchmark USPTO-50k. Next, we investigate how MEGAN scales to the large-scale retrosynthesis task. Finally, we present results on forward synthesis prediction on the USPTO-MIT dataset.

Retrosynthesis
Data First, we evaluate MEGAN performance in retrosynthesis prediction. The goal is to predict correctly the set of reactants based on the product of a reaction. The accuracy is measured by comparing the SMILES 18 representation of the predicted molecules to the SMILES representation of the ground truth target molecules. Before the comparison, we remove atom mapping information from the SMILES strings and canonicalize them using RdKit 45 . We use top-k accuracy computed on reactions from the test set as the main evaluation metric.
We evaluate on the USPTO-50k dataset of approximately 50000 reactions, which was collected by Lowe 49 and classified into 10 reaction types by Schneider et al. 50 . We use the same processed version of the dataset as Dai et al. 13 , where each reaction consists of a single product molecule and a set of one or more reactants, with corresponding atoms between the reactants and the product mapped. Following other studies, we split the reactions into training/validation/test sets with respective sizes of 80%/10%/10%. We use the dataset and split provided by Dai et al. 13 .  51 . We excluded from the comparison some of the concurrent works which had an information leak problem in evaluation (see: Section Concurrent methods excluded from the comparisons in the Supplement).
Results Table 2 reports results on the USPTO-50k benchmark in a variant with and without reaction type information provided. When the reaction type is given, MEGAN outperforms all baselines in all metrics, with the exception of top-1. When the reaction type is unknown, MEGAN beats prior models for K ≥ 3 (and concurrent models for K ≥ 10) and achieves comparable results for lower values of K.
We hypothesize that the advantage of MEGAN for large K stems largely from the fact that MEGAN generates reaction as a sequence of edits. This might help to efficiently search through different plausible reaction centers, hence covering a more diverse subset of the reaction space. It can also enable MEGAN to achieve high coverage of the reaction space, which is indicated by the top-50 accuracy of 93.2% when reaction type is unknown and 95.3% when the reaction type is provided.

Large scale retrosynthesis prediction task
Data The USPTO-50k benchmark, although of relatively high quality, is a small dataset containing only 10 specific types of reactions. We train MEGAN for retrosynthesis on a large benchmark to test its scalability. We use the original data set collected by Lowe 49  Results Table 3 shows top-k accuracy on USPTO-FULL for MEGAN compared to other methods for retrosynthesis prediction. We see that our model achieves competitive performance on a large scale retrosynthesis data set, slightly outperforming other methods in terms of top-10 accuracy. identifies pairwise atom interactions in the reaction center and ranks enumerated feasible bond configurations between these atoms. WLDN5 31 improves this method by combining the problems of reaction center prediction and candidate ranking into a single task. GTPN 23 predicts actions on the graph of substrates, similarly to MEGAN, however, it limits them to actions between existing atom pairs. We include also the concurrent AT model 51 in the comparison.
Experimental Setting We remove the CHIRAL TAG and BOND STEREO features from the model input, as USPTO-MIT has no stereochemical information. Forward synthesis usually takes fewer modifications to predict than retrosynthesis, so we reduce the maximum number of steps to 8 and use a beam size of 20 for evaluation. We also change the strategy to break ties in action ordering. For forward synthesis, we prioritize actions that act on most recently modified atoms. Apart from these changes, we use the same architecture, hyperparameters, and training procedure as for USPTO-50k.
Similarly to Schwaller et al. 15 , we train for two variants of forward synthesis prediction. For the separated variant, compounds that directly contribute to the product are explicitly marked in the set of substrates with an additional atom feature. For the mixed variant, such information is not provided, so the model has a harder task as it has to determine the reaction center from a larger number of possible reactants.
In forward synthesis, MEGAN can generate multiple products for a given set of input substrates set, whereas the USPTO-MIT set contains single-product reactions. We use a simple heuristic to acquire the main reaction product from the MEGAN prediction by selecting the product that has the longest SMILES string. From our observation, this heuristic works correctly for most of the reactions and does not lower the Top K accuracies.
Results training strategy for forward prediction is a promising topic for future work.
In Table 5 we compare the performance of MEGAN and Molecular Transformer on USPTO-MIT in top K predictions. We use the best non-ensemble models provided by Schwaller et al. 15 and set the beam size to 20 during evaluation. We observe that MEGAN surpasses the accuracy of Molecular Transformer for high K values, which again indicates its ability to explore the reaction space efficiently.

Analysis and ablation studies Error analysis
To better understand model performance, we perform here an error analysis focusing on retrosynthesis prediction on the USPTO-50k dataset.
We first visualize in Figure 2 MEGAN prediction on the USPTO-50k test set for random 5 target molecules on which the predicted substrates are different from the ground truth reaction. We observe that the first four reactions are feasible and can be executed using standard methods. The last reaction may pose difficulty because of the nontrivial regioselectivity that depends on the exact reaction conditions utilized.
To better quantify the phenomenon that the top proposed reaction is often correct despite being different from the ground truth, we sampled random 100 reactions from the USPTO-50k test set for which the predicted substrates by MEGAN differed from the ground truth substrates. We then analyzed the correctness, by labeling them into the following categories: (1) no issues detected (seen as a true chemical reaction by a chemist), (2) incorrect chirality of the substrates, (3) the reaction has a low yield and/or important side products, (4) the reaction can be executed but only in multiple stages (e.g. with use of protecting groups), (5) there is a reactive functional group that participates in the reaction instead, (6) incorrect reaction for another reason. All 200 reactions were shuffled and assigned randomly to two experienced organic chemists.
We show the break-down of MEGAN errors by these categories in Figure 3 and examples of reaction in each category in Figure 3. On the whole, in 79.6% cases, the top prediction by MEGAN was deemed correct by chemists even though it differed from the ground truth prediction.
In comparison, 89.5% of ground truth reactions were considered correct. Incorrect chirality of the substrates was the most common reason for reaction to be labeled as incorrect (7.1%, compared to 0.5% of ground truth reactions).
The main outcome of this error analysis is that for retrosynthesis prediction on the USPTO-50k, MEGAN performance seems to be close to human-level performance in the sense that the proposed reactions start to be hard to distinguish from ground-truth reactions in terms of their correctness.
Having said that, there is still significant room for improvement, for example in terms of the chirality of the proposed substrates. Another common source of errors was ignoring the existence of a reactive functional group in the substrates.

Ablation study on action ordering
Our gradient-based method of training MEGAN is based on a fixed ordering on ground-truth actions. Defining this ordering requires breaking ties between different actions of the same type.
For example, if there are several EditBond actions, we have to decide which should be executed first by the model.
In retrosynthesis prediction tasks, we broke ties based on the last time a given atom was modified.
Actions that modify atoms that were least frequently modified are prioritized. The remaining ties were broken at random. We refer here to this ordering algorithm as BFS RAND-AT.
We explore here two modifications to this ordering. First, instead of prioritizing actions that Figure 2: The top-ranked retrosynthesis prediction by MEGAN for 5 random target molecules from the USPTO-50k test set on which the prediction is different from the ground truth reaction. The first four reactions are feasible and can be executed using standard methods. The last reaction may pose difficulty because of the nontrivial regioselectivity that depends on the exact reaction conditions utilized.   Figure 3: Categorization of the top-ranked retrosynthesis predictions by MEGAN on a subset of USPTO-50k test set on which the top prediction differs from the ground truth reaction. For reference, we also include labels assigned to ground truth reaction from the test set. We find that the top-ranked prediction is often (in 79.6% cases) deemed correct by a chemist even though it differs from the ground truth ("No issues detected"). The most common source of error is incorrectly predicting chirality in the substrates ("Incorrect chirality", 7.1%).  Figure 3).

Target product
Ground truth reactants

Predicted reactants
Only possible in multiple steps Low yield or side products Incorrect chirality A reactive functional group ignored modify the least frequently modified atom, we experiment with prioritizing actions that modify the most frequently modified atom. We refer to this as DFS ordering. Second, we investigate a strategy to break remaining ties based on the canonical SMILES. In this strategy, we prioritize actions that modify atoms that are earlier in the ordering determined by the canonical SMILES. We refer to this as CANO-AT. By combining these different choices, we arrive at four different ordering strategies: BFS RAND-AT, DFS RAND-AT, BFS CANO-AT, and DFS CANO-AT.
We present the performance of all variants on the USPTO-50k dataset in Table 6. We observe that BFS RAND-AT ordering achieves the highest performance across most values of K, which motivates our choice to use it for training MEGAN on retrosynthesis prediction tasks. It is worth emphasizing that the choice of action ordering can strongly impact the performance of MEGAN.
For example, the RANDOM action ordering achieves over 4% lower top-1 accuracy compared to DFS RAND-AT, which we used in the experiments. Performance with respect to reaction popularity Next, we analyze how MEGAN performance depends on reaction popularity. To approximate reaction popularity, we count the number of times that the corresponding template occurred in the training set. We use code provided by Coley et al. 31 to extract templates for each ground-truth reaction from the datasets. Crucially, MEGAN generates reaction as a sequence of graph edits, which is inspired by how chemists represent chemical reactions. This inductive bias might explain the strong empirical performance of the model. In particular, we argued that it enables the model to more efficiently search through the space of plausible reactions. Looking forward, generating reactions as a sequence of edits is promising for building more intuitive human-computer interaction in synthesis planning software.
A natural topic for the future is reducing the reliance of MEGAN on the mapping between products and substrates. We also found that incorrectly predicting chirality of the output was the most common source of errors on the USPTO-50k dataset. Improving these aspects has the potential to further push state-of-the-art in retro and forward synthesis prediction.

Data and Software Availability
We open-sourced our code and trained models at https://github.com/molecule-one/megan. Our implementation is based on freely available software and is released under the MIT license. All datasets are available online, with references given in the main text.

Supplementary Material Featurization
We featurize atoms and bonds with one-hot encoded vectors of features calculated using RdKit 45 .
On each dataset, we select features that allow for correct reconstruction of all products and substrates from the development (training+validation) set. In Tables   For both atoms and bonds, we add a special IS EDITED feature that marks all bonds and atoms that have been modified by actions. This aims to help the decoder to focus on these atoms and bonds, as they are the most probable candidates for the next generation steps.

Graph edit actions
In Table 8 we show all possible graph actions found on the USPTO-50k development set. The actions were found during the generation of the training samples, which we describe in the next section. The actions have different sets of parameters, depending on the action type. BondEdit actions are bond actions, that is they are predicted for a pair of atoms. All other types of actions are atom actions and are predicted for a single atom in the graph.

Ground truth action ordering
An important design choice in training of MEGAN is the algorithm that determines ordering of ground-truth actions that transform the product into the substrates (or vice versa in forward prediction). This ordering allows us to backpropagate through the maximum likelihood objective on the training set.
Action type priorities At first, we define a general ordering of action types for forward-and retrosynthesis in Table 1. This ordering defines which types of actions should be performed first if there are more than one type of actions that could be applied to an atom in a generation step. For retrosynthesis, we give bond deletion the highest priority, as it is usually the step that determines the reaction center, which is the first step of reaction prediction in other methods. [36][37][38]48 Analogously, adding a bond has the highest priority for forward synthesis, as this usually determines the reaction center in forward prediction. Priorities of other types of actions were determined experimentally; we did not observe a significant change in validation error when modifying these priorities.
Atom priorities Given the ordering of action types described in Table 1, we still need to decide on which action to choose in a scenario when two of the ground-truth graph edits are possible for different atoms. For example, we might need to decide on the order of adding leaving groups to two atoms that were disconnected by deleting a bond between them. Ordering of such actions can be determined by two factors: 1. Whether an atom has already been modified (BFS vs DFS) 2. Order of the atom in the canonical SMILES of the target compounds (CANO-AT vs RAND-AT) In the case of EditBond actions that act on pairs of atoms, their priority is determined based on the highest of the priorities of these atoms. We present the whole algorithm for determining the next action in the Algorithm 1. Let Act(a) = all actions that lead to T that could be performed on a Let Act (a) = all actions from Act(a) with the action type that has the highest priority Let ActionT ype = action type of actions from Act (a) if ActionT ype = EditAtom then By definition, there could only be a single EditAtom action performed on an atom. return The single EditAtom action from Act (a) else if ActionT ype = EditBond then Let b = the second atom from actions in Act (a) ranked highest according to the atom ordering return The EditBond action between a and b else if ActionT ype = AddAtom then Let b = the second atom from actions in Act (a) ranked highest according to the atom ordering return The AddAtom action between a and b else if ActionT ype = AddBenzene then return A random AddBenzene action from Act (a) end if Hyperparameter search Table 9 shows the hyperparameter values used for the final model on the USPTO-50k. We found these hyperparameters manually based on the performance on the USPTO-50k validation set. For other datasets, we only slightly modified these hyperparameters (as described in the main text).

Concurrent methods excluded from the comparisons
After a detailed examination, we decided to exclude from the benchmark comparisons some of the concurrent methods of retrosynthesis prediction that suffer from an information leak problem between the training and test sets. This problem is described in details at https://github.com/utasmile/RetroXpert (accessed on May 7th, 2021). The methods with the information leak problem include GraphRetro 37 and RetroXpert. 38