extfieldsofstudy:
Physics
Computer Science

text:
QuISP: a Quantum Internet Simulation Package

We present an event-driven simulation package called QuISP for large-scale quantum networks built on top of the OMNeT++ discrete event simulation framework. Although the behavior of quantum networking devices have been revealed by recent research, it is still an open question how they will work in networks of a practical size. QuISP is designed to simulate large-scale quantum networks to investigate their behavior under realistic, noisy and heterogeneous configurations. The protocol architecture we propose enables studies of different choices for error management and other key decisions. Our confidence in the simulator is supported by comparing its output to analytic results for a small network. A key reason for simulation is to look for emergent behavior when large numbers of individually characterized devices are combined. QuISP can handle thousands of qubits in dozens of nodes on a laptop computer, preparing for full Quantum Internet simulation. This simulator promotes the development of protocols for larger and more complex quantum networks.

Introduction
The second quantum revolution has steadily been gaining momentum over the last two decades [19]. Quantum computers are at the forefront with private companies racing to build ever larger quantum devices demonstrating quantum supremacy [3,68]. Quantum networks [58] promise to revolutionize the field of communication by exploiting the rules of quantum mechanics in order to bring enhanced, and in some cases completely new, functionality. Some of the main promises of quantum networks include distribution of secret keys for secure communication where a malicious party is doomed to be discovered due to the fundamental laws of quantum mechanics [6,23]. The possibility of distributed quantum computation where quantum computers are networked together in order to solve a difficult computational task in a cooperative fashion [15] is being explored. Blind quantum computation [9] allows a client with limited quantum resources to delegate their computation to a powerful quantum server without revealing the input, the computation itself or its output. Enhanced clock synchronization [29,35] promises to improve the accuracy of global navigation while distributed quantum states can be used to construct ultra-sensitive sensor networks [53,69]. These applications are likely just the tip of the iceberg and novel uses of quantum networks will be discovered.
The fundamental resource behind these marvelous applications is entanglement [27], the ability of spatially separated quantum states to be correlated more strongly than classical states. The primary job of a quantum network is to distribute these entangled states between two or more parties. The ultimate goal is to build a global Quantum Internet [36,56,64] of entangled quantum devices, all running on heterogeneous hardware yet still being able to engage in reliable and efficient quantum communication.
Quantum networks are quickly becoming more than just a dream, with early quantum key distribution (QKD) networks having been implemented already [11,22,24,34,51,55,57]. Currently, the race is on to demonstrate feasibility of repeaterbased quantum networks, where entanglement over large distances is constructed recursively from shorter entangled links [8,20], with some of the basic hardware elements being demonstrated experimentally [4,32,52].
Besides the immediate technological challenges facing us when developing a global quantum network, we have to also deal with open questions for the future of the Quantum Internet. This is where a quantum network simulator becomes a crucial research and design tool. Some areas where a good simulator will be invaluable are: (1) Protocol design. Particularly testing of detailed protocol design to validate correct operation and study of the interaction between classical and quantum portions of the network. (2) Connection architecture and performance prediction. Proposed generations of quantum networks [44] exhibit complex behavior making analytic prediction of their performance difficult with realistic parameters. (3) Dynamic behavior. Important open questions are the stability of quantum networks as conditions change over time and the response of protocols to dynamically changing network topology. (4) Emergent behavior. As the system increases in scale, new and unexpected behavior may emerge that current naive models of quantum networks do not take into account. How do competing connections behave? In particular, are quantum networks subject to congestion collapse, or can short-distance connections starve long connections? Must we trade connection fidelity (quality) for performance? Naive models suggest that end-to-end, high fidelity connections are always possible, but it is still an open question whether this is true in a dynamic, global network.
Our approach to simulating large-scale quantum networks using our Quantum Internet Simulator builds on OMNeT++ [63], a modular, component-based architecture simulation environment. Its focus is on protocol design for complex, heterogeneous networks at large scale while keeping the physical layer as realistic as possible. Following on from earlier simulations of single lines of repeaters [30,59] and an early network simulator [2], the past two years have seen a flurry of activity in the field of quantum network simulation, with the introduction of a number of simulators such as NetSquid [14], SeQUeNCe [66] and QuNetSim [18]. These simulators are limited by focusing on physically realistic simulation of a single, small network. QuISP is designed with internetworking in mind while maintaining full physical realism. Our long-term goal for the simulator is to be able to handle an internetwork with 100 networks of 100 nodes each, with each network running independent error management protocols, hardware parameters, and topology.
The manuscript is structured as follows. First, we provide basic concepts of quantum information processing and our new quantum state representation proposal towards to large quantum simulation in Section 2. We then introduce the background design of quantum network architecture in Section 3. Our main proposal for our quantum networking simulator is described in Section 4. We give the detailed explanation of basic design principles and the current implementation of QuISP. Finally, we demonstrate experiments to show the correctness and performance of QuISP in Section 5 before concluding in Section 6. Additional details on configuring the simulator and the mathematics of our new error basis simulation are in the appendices.

Quantum
In this section, we give a brief overview of the basics of quantum information processing and discuss our approach to scalable simulation of quantum networks.

Quantum Information Processing
Quantum networks encode information into quantum states [46]. The simplest and also the most commonly used quantum system is a qubit. It is comprised of two basis states, a |0 (pronounced "ket 0") and a |1 (pronounced "ket 1"). Unlike a classical bit that is either a 0 or a 1, the qubit can be in a superposition of the two base states, written as α|0 + β|1 , where α and β are known as complex probability amplitudes. They give the relative likelihood that the state of a qubit is a |0 or |1 . In particular, the state |+ = (|0 + |1 )/ √ 2 is an equal superposition of |0 and |1 . This superposition principle extends to multiple qubits as well. In particular, qubits A and B can be in a superposition α|0 A |0 B + β|1 A |1 B . This is a rather special state because it is impossible to write it as a product of local qubit states of A and B. Such states are known as entangled states and are at the heart of most quantum technologies. A particular entangled state used ubiquitously in quantum communication is when α = β = 1/ √ 2, known as a Bell pair.
The basic way to transform the state of a qubit is via Pauli operators X, Y , and Z. Pauli X is the quantum analogue of a bit-flip. It flips the state of a qubit X|0 = |1 and vice versa. Pauli Z does not have a classical counterpart as it introduces a phase, namely Z|1 = e iπ |1 = −|1 . This phase is important when the qubit is in a superposition state since Z|+ = (|0 − |1 )/ √ 2 is experimentally distinguishable from |+ . Quantum state tomography [1] is the process of characterizing the quantum state by measuring identical copies in different bases in order to reconstruct a statistical description of the state known as the density matrix. It is key to extracting useful information about the quantum state, particularly its fidelity which is a measure of how close the actual output of a quantum protocol is to its desired one.
The state of a qubit can be transmitted in a quantum network using quantum teleportation without physically sending the system encoding the qubit [7]. Consider a client wishing to communicate the state of qubit CQ 1 to a server. They must share two qubits in a Bell pair, denoted by CQ 2 and SQ. The client performs a Bell-state measurement (BSM) on qubits CQ 1 and CQ 2 , and sends the classical result of this measurement to the server who performs a conditional correction on its qubit SQ. This ensures that the state of qubit SQ is the same as that of qubit CQ 1 was initially. The entanglement between qubits CQ 2 and SQ is consumed during the measurement and therefore must be reestablished prior to any further quantum teleportations. Note that entanglement does not imply any directionality so the server can transmit its qubit back to the client provided they share a Bell pair.
Quantum teleportation can be extended to pairs of entangled qubits, in which case it is known as entanglement swapping [28]. Consider two entangled pairs, denoted by A-B 1 and B 2 -C as depicted in Figure 1. Measurement in the Bell basis of qubits B 1 and B 2 will result in qubits A and C becoming entangled. Entanglement swapping is the basic building block of quantum repeater networks as it allows end-to-end entangled connection between distant nodes to be established.
Sharing entanglement between neighboring nodes of a quantum network is not a deterministic process due to noise Dirac's bra-ket notation sparse representation of state vector (also used for variable names) O(2 n ) but grows sharply to O(2 n ) von Neumann's density matrix pure or mixed (error) states (sparse representations also possible) O(4 n ) Stabilizer [26] shorthand for specific states, specifies constraints on states (eigenoperators) Tensor network [40,48]   and losses in the fiber. Even once the entanglement is established its quality deteriorates over time due to decoherence of quantum memories. One of the jobs of a quantum network is to manage these errors via entanglement purification [5]. In this protocol, two imperfect entangled pairs can be converted probabilistically into a single entangled pair of higher fidelity.

Simulating in the Error Basis
Quantum systems are exceptionally fragile and susceptible to errors due to interaction with the environment. Furthermore, distant nodes in a quantum network communicate via exchange of single photons, meaning photon loss in the optical fiber becomes a major source of errors. It is therefore crucial that a quantum network simulator models all physically relevant sources of errors accurately.
Calculating the state of a quantum system by hand or classical computer can be done using one of several representa-tions, with differing purposes and scaling, as shown in Table 1. Many important states can be written down using shorthand methods such as the sparse Dirac bra-ket notation or stabilizers [26], but complete representation of an arbitrary state may require up to 2 n complex numbers when the state is pure, or error free. Using special techniques on classical supercomputers, computations of special cases involving up to 56 [50] or 81 [10] qubits have been simulated.
Traditional tracking of the evolution of a mixed, or noisy, quantum state composed of n qubits is substantially harder. The density matrix can represent any quantum state after any computation or noise process, but when fully populated, requires 2 n × 2 n complex numbers. In the context of quantum networks, full density matrix simulation is tractable for 1G networks, where the simulation must track many independent quantum states, but each state is composed of at most four qubits. This approach has been used in numerous simulations, especially those focused on tuning low-level hardware control parameters [39,49]. However, this approach quickly becomes intractable when moving to 2G and 3G networks, where logical Bell pairs are encoded into a large multi-qubit state for quantum error correction [44]. Simulating only two 2G links using a 7-qubit error correcting code on each link [31] requires 16 × 2 56 bytes, about one exabyte. Simulating a 100 hop long connection could take a mind-boggling 3.2 × 10 3011 bytes of memory. Thus, it is imperative to separate simulators based on their purpose; our goal is to study the architecture and protocols of networks, rather than the physics of devices or complete quantum algorithms.
QuISP works on the premise that the desired quantum state is known and we only have to track deviations from this ideal state in the form of errors that are affecting it. This approach is adapted from quantum error correction [17,45] and represents a novel way of simulating quantum networks. Deviations from the ideal state can be tracked efficiently leading to scalable simulation of truly global quantum networks beyond 1G for the first time. This scalability is achieved by exploiting the fact that quantum networks only need to apply a fixed number of operations and the major sources of errors can be discretized into a small set. Our approach is tailor-made for quantum network simulation and is not suitable for simulation of a universal quantum computer.
The simulator tracks Pauli X, Y , and Z errors discussed in Section 2. It also tracks relaxation and excitation errors, denoted by R and E, respectively. Relaxation captures the process of energy loss due to environmental noise by transforming any initial state to |0 , while excitation captures the opposite process where any arbitrary initial state is transformed to |1 . Finally, the simulator tracks photon loss L as well. This list does not form an exhaustive set of all possible possible errors and may be expanded or shrunk depending on the particular physical scenario under consideration.
The quantum state at time t is represented by a m + 1element error probability vector π(t), where m is the number of errors that we are tracking. Each element π j represents the probability that the quantum state has been affected by a particular error j at time t, The error probability vector is normalized, meaning ∑ j π j = 1. The evolution of the error probability vector is given by the transition rate matrix Q. It is a right stochastic matrix with elements Q i j representing the probability per unit time of transitioning from error state π i to error state π j , and they satisfy the normalization condition ∑ j Q i j = 1. The error probability vector π(t) a time t is given by The transition rate matrix Q is independent of time and the error probability vector π(t) depends only on π(t − 1), meaning the evolution of the qubit is modelled as a Markov process. Detailed discussion of the transition matrix can be found in Appendix B. The simulator samples the error probability vector π(t) prior to an operation on the qubit and applies the sampled error type 1 . The simulation must be repeated to gather statistics about the real state of the qubit and gain information about the qubit's fidelity.
Extension to multi-qubit systems is straightforward. The state of N qubits in QuISP is described by an N(m + 1)element error probability vector where the first m + 1 entries describe qubit 1, second group of m + 1 describe qubit 2 and so on. The full transition rate matrix for N qubits is given by a block-diagonal matrix composed of N single-qubit transition rate matrices. 1 The current version of the simulator samples the error probability vector before a measurement or purification is performed. The next version will extend the sampling to entanglement swapping as well.

Network Design
In this section, we discuss the basic network design that QuISP assumes.

Quantum Network Architectures
There is currently no consensus on the best overarching network architecture for quantum networking, though the key principles are coming into view [38,61] and some protocol elements have been proposed [16,37,41]. Supporting further research in this area is the primary purpose of QuISP. However, it is becoming clear that some basic hardware and software components will most likely be shared between future candidate architectures. We give a brief outline of these components in this subsection.
Hardware architecture: There exist a number of potential candidate physical systems that are suitable for encoding qubits, broadly divided into two categories. Stationary qubits or matter qubits are envisioned to store and process information at the nodes of a quantum network, acting as the hosts. Candidate physical systems include nitrogen-vacancy centers in diamond [52], trapped ions [21], atomic ensembles [67] and superconducting qubits [43].
Inter-node quantum communication is achieved by using flying qubits encoded onto single photons travelling through optical fibers. We refer to these fibers as quantum links. Photons are ideal information carriers as they do not interact strongly with their environment and they travel at very high speeds. Using flying qubits, it is possible to distribute entangled Bell pairs between two distant nodes of a quantum network. This can be achieved using one of the three existing quantum link architectures [33] depicted in Fig. 2. Memory-Memory (MM) link connects two quantum nodes directly where either node is equipped with a Bell-state analyzer (BSA), an optical device that performs a Bell-state measurement on two incoming photons. Memory-Interfere-Memory (MIM) link places the BSA in the middle of the quantum link. Memory-Source-Memory (MSM) link replaces the the BSA in the middle of the quantum link with a source of entangled photonic pair states (EPPS). Despite the architectures appearing fairly similar, they differ significantly in their performance as well as the technological maturity required to implement them.
Since quantum communication operates at the singlephoton level, attenuation becomes a major source of error. Unlike classical bits, qubits cannot be copied and resent owing to the no-cloning theorem [65], a fundamental property of quantum mechanics forbidding to deterministically copy arbitrary states of quantum systems. Amplification at the level of single photons becomes ineffective as well [12,13]. This limits the practical length of quantum links to mere tens of kilometers.
In order to get around this problem, a new type of node was introduce known as a quantum repeater [8,20]. One of the jobs of a quantum repeater is to share Bell pairs with its neighboring nodes and implement entanglement swapping in order to create a Bell pair between these nodes. In this way, the no-cloning theorem can be sidestepped and photon loss mitigated, resulting in the possibility of establishing end-toend Bell pairs between arbitrarily separated quantum hosts. Figure 3 depicts the hardware and software components of a quantum repeater.
A particularly important component of the quantum repeater is the Quantum Network Interface Card (QNIC). The QNIC is the quantum analogue of a classical NIC with one major difference being that a QNIC is able to apply quantum operations to the store quantum information, making it a quantum computer with limited capabilities. In particular, the QNIC is capable of applying single-and two-qubit gates as well as single-and two-qubit measurements.
Software architecture Classical software running on a quantum repeater will play a crucial role when designing efficient repeater-based quantum networks. Our proposed software architecture, Quantum Repeater Software Architecture (QRSA), as shown in Figure 3, consists of five software components.
•  These components communicate as needed to convey when to start operations and the current status of devices.
An end node has almost the same functionality as a repeater, but it also has an Application component responsible for performing end-to-end applications.

RuleSet Protocol
In order to distribute end-to-end entanglement, both end nodes and quantum repeaters must know what actions to perform, when to execute them, and what other nodes are taking part in the process if the actions need to be coordinated. For example, the repeater must know the nodes that it shares Bell pairs with when executing entanglement swapping since the results of the procedure must be shared with those nodes.
To this end, the RuleSet protocol was proposed in [41], which QuISP supports. The goal of this protocol is decentralized, autonomous but coordinated actions of the quantum repeaters with minimal classical inter-node communication. Figure 4 is an example of RuleSet structure. The RuleSet is a collection of Rules such purification and entanglement swapping. These RuleSets are built in the connection setup phase discussed in Section 3.3, and executed in a specified order. After being acted upon, the Bell pairs belonging to a particular Rule are passed to the next in sequence.

Rule 2
Resource Table   Action Resource sources (Bell pairs). In quantum networking, shared Bell pairs must be managed by each node in a coordinated fashion and appropriately structured RuleSets provide this required consistency in terms of quantum operations.
Condition Clauses are composed of single or multiple conditions to be met before the Action is executed. For example, if node A requires two entangled states with node B to perform one action, A must track the number of shared entangled states with B. In such a situation, the Condition used is the Enough Resource Clause, which is satisfied when the number of total entangled pairs shared with the proper partner is larger than a threshold (In this case, the threshold is two). Other than Enough Resource Clause, there are several clauses supported in this simulator.
An Action Clause is a set of operations including resource assignment changes, qubit manipulation, and classical message transfer. Once Condition Clauses are met, the corresponding action is immediately executed. For example, Swap refers to the resource table that belonging to the Rule and recognizes the corresponding qubits. Then, this action chooses one state entangled to its left and one entangled to its right and applies Bell state measurement. Informing the partners of the result of the Bell state measurement is one responsibility of the action.

Connection Setup Request
Accumulate interface information of intermediate nodes

Initiator Responder
Requirements Interface info Link cost Figure 6: Connection setup process to establish agreements between initiator and responder

Connection Setup
The connection setup is the step requires to gather all the information to create RuleSets to be executed for nodes that will be participating in the end-to-end Bell pair generation.
The connection setup process used in QuISP is adapted from protocol outlined by Van Meter and Matsuo [60]. Figure 6 shows the procedure of the connection setup. It involves a two pass process, gathering the link information along the path starting from the node that tries to establish the connection (Initiator) and planning the RuleSet to be distributed among the nodes along the path at the other half of the connection (Responder). The first part, at the Initiator node, it receives the requirement for the connection from application level, like the quality of the connection (fidelity of the end-to-end Bell pair) and the number of Bell pairs. In this outbound pass, every node along the path will include their link characteristics into the message, reserve the QNIC, and relay this connection setup message to the next hop. If the QNIC cannot be reserved because it is already in use for another connection, the node will reject the request and the connection setup reject message will be sent to all the previous nodes along the path.
When the connection setup message arrives at the Responder node, the Responder's job is to plan out how each node should execute their share of work in order for the end-to-end Bell pair creation to succeed. After planning out and creating RuleSets for all nodes, the Responder sends Connection Setup Response with the RuleSets back to all nodes along the path.

The Simulator
In the previous section, we explained the network design underlying QuISP. Here we introduce its implementation and use, with an example simulation.

Goals and requirements
The long-term goals for developing QuISP are supporting 1G, 2G, and 3G quantum networking (generations classified by the type of error handling) and internetworking in a scalable manner to establish highly reliable protocols for each generation. In the higher generations of quantum networking, thousands of qubits will work together when creating a single end-to-end logical Bell pair. Furthermore, this simulator will allow us to replicate the dynamic behavior of quantum networks and internetworks with complex network topologies, including the concepts of heterogeneity and network boundaries.
As noted in the introduction, our long-term goal is 100 networks of 100 nodes each. With inter-repeater spacing as low as 10km for some hardware architectures, long chains of repeaters between routers will be necessary to achieve distance. With low initial link performance, the router-level topology will have to be rich in order to provide adequate connectivity and resource availability. We expect this scale of simulation to aid in assessing an internetwork with such real-world constraints.

Basic Design Principles
We now highlight the main principles that we used as a guide in designing QuISP.
Realism. The simulation should provide accurate information about the physical states that are being distributed in the quantum network. This in turn provides accurate information about the fidelity of the distributed states. Emphasis in the design of QuISP has been placed on realistic noise models. Qubits stored in quantum memories will undergo decoherence processes that wash away their quantum properties and result in deteriorating fidelity. Photon loss in fiber is a major source of error and therefore must be accounted for in the simulator. Early simulation work assumed 2 n hops, all of the same length and fidelity [30,59], but that is not realistic.
The implementation of the physical layer simulation must be kept cleanly separated from the implementation of software for the nodes themselves. The router software can know only what it would know in the real world: information it can learn from classical measurement outputs from the quantum hardware and from exchanging messages with other nodes.
Scalability. One of the main goals of this simulator is the study of large-scale quantum networks and their emergent complex behavior. Scalability concerns come in two flavors: the number of qubits involved in a single quantum state, and the number of nodes, links and connections in the whole network. General multi-qubit states are difficult to simulate classically. Our novel approach of working in the error basis discussed in Section 2.2 allows us to simulate quantum networks at the cost of only a few classical bits per qubit, and makes QuISP suitable for simulation of truly large-scale quantum networks of any generation.
Flexibility. This simulator has been designed to offer the user a great deal of customizability. Currently, there exist a number of physical systems that are considered to be promising candidates for implementing the physical layer of a quantum network. However, all these candidate systems differ substantially in terms of parameters such as memory coherence times and operation clock rates. QuISP allows the user to customize the error types that affect the physical layers as well as the individual parameters corresponding to error rates. This ensures that QuISP is capable of accurately simulating any current and future quantum hardware. Furthermore, the topology of the quantum network can be defined easily and this simulator comes with a set of predefined number of qubits per node and topologies as well.

Event Simulation with OMNeT++
OMNeT++ is a C++-based discrete event simulator primarily for classical networking simulations [47,54,62,63]. The OMNeT++ model consists of modules that communicate via messages. The first batch of messages is created during module initialization. Their destination can be any other modules connected to the sender module or to self. Upon receiving a message, which constitutes an event, the module processes the message and usually creates another batch of messages to be sent. The only proper way to trigger execution of a function in a certain module is by receiving messages. It is assumed that the processing of events takes zero time and that nothing of significance occurs between two events. The simulation ends when there are no messages left to be processed or the simulation reaches the time limit set by the user. The simulation time is tracked by the timestamp of events, which is determined by the time a message is scheduled to be sent.
Modules and messages in OMNeT++ are provided by the cSimpleModule and cMessage class, accordingly. cSimpleModule is the active component in OMNeT++ and its functionality can be extended by adding C++ code or by composing multiple simple modules into a compound module.
Furthermore, QuISP also provides a graphical user interface from the OMNeT++ IDE to manipulate the network topology stored in a NED file, which defines network devices in very structured manner. An example of a NED file is shown in Appendix A. The simulator also allows the user to change the initial simulation parameters or variables on launch using an INI configuration file. It is also possible to visualize many stats, such as the changes of parameters during the simulation. These functions allow us to debug and check if the simulation works properly and to make intuitive performance measurements.

Basic Implementation
A concise diagram outlining our implementation is shown in Figure 7. Software modules (QRSA, Application, etc.) are built on top of cSimpleModule, which supports standard message handlers. Stationary Qubit is also a subclass of cSimpleModule, but it holds its quantum properties as its variables and behaves as quantum memory. Those software components can emit registered classical packets such as ConnectionSetupRequest written in cMessage. Flying qubits can be thought of as an emitted message from Stationary Qubit. Thus a flying qubit is implemented as an extension of a classical packet passed through quantum channels.
A quantum channel as a carrier of the flying qubit is implemented as a subclass of cDatarateChannel. The errors affecting flying qubits are calculated based on the distance between two nodes.

Supported Functionality
QuISP currently supports end-to-end tomography-based quality analysis of connections as the application. The output of the simulation is the characteristics of the path used for each connection, explained in detail in Section 5.
The user can define complex topologies easily. Almost every factor that plays a role in the link characteristics such as node placement, distance between nodes, link types, number of buffer qubits in each interface, channel error and loss rate, and the efficiency of most components that contribute to 1G and 2G technologies can be set individually or from defaults. The simulator also supports options for selecting how the end-to-end Bell pair is generated. The order of entanglement swapping can be chosen between the usual binary tree like structure [8] or swap-once-ready at each node. The user can also tinker with the number of rounds of purification performed before entanglement swapping takes place and choose one of the four types of purification methods [8,25] we currently provide.
Since most of the groundwork for the RuleSet-based protocol has been implemented, curious users can modify the simulator code, write their own entanglement swapping policy or purification method, and experiment with how to achieve higher quality connections.

Configurable Parameters
As we mentioned in the previous section, users can make their decisions on the error rate, the number of purifications, the way of entanglement swapping. Currently, there are more than 50 types of parameters that can be given through the INI file. The following parameters are most likely to be configured by users.
• Network. The name of the network written in the NED language, specifying the number and types of nodes and the network topology.
• Channel errors. The errors on a quantum channel. This specifies how much and what types of errors in a particular distance. Currently, all Pauli errors and photon loss errors are supported, and users can specify the ratio of each error.
• Memory errors. In addition to Pauli errors, we support relaxation and excitation errors. As the simulation time passes, this error accumulates on each qubit.
• Gate errors. In a realistic situation where we assume that devices are imperfect, there are always small operational errors in quantum gates. We can handle errors on basic gate sets (Pauli gates, Hadamard, Controlled-X).
• Measurement errors. Quantum memories must be measured to get internal information which affects the quantum state itself. Measurement errors are composed of Pauli errors.
• BSA errors. It is possible to have different errors for internal and external BSAs. Tunable parameters are photon loss, photon detection rate, and dark count probability.
• Photon emission success probability. The probability that a quantum memory successfully emits a photon which is then captured by the fiber.
• Traffic pattern. A configurable parameter for traffic generation pattern. Currently, this simulator supports single traffic generation from one node to the other random node and multiple traffic generation from all end nodes to randomly selected nodes.
• Measurement count. The number of measurement in tomography. When the required precision of tomography is higher, the number of measurement needs to be larger.

Simulation Output
The simulator returns the performance statistics at the end of simulation when the target configuration contains tomography. The output files contain the following information.
• Fidelity. Calculated by the software tomography process, this indicates the quality of Bell pairs between two end nodes.
• Bellpair per sec. The number of Bell pairs generated in one second (in simulation time).
• Tomography time. The time until the entire tomography process finishes.
• Tomography measurements. The number of measurements of tomography.
Besides these properties, the output statistics include more link information. An example output is given in Listing 2 in Section 4.8.

Example Simulation
We demonstrate the workflow to simulate end-to-end entanglement generation between Initiator and Responder on the network shown in Figure 8. First, we prepare a NED file corresponding to the target network. (See Appendix A for more details.) The nodes are defined as submodules of the network. After enumerating all involved nodes, we define both classical and quantum connections between these nodes.

Listing 1: Example configurations in INI file
Second, we prepare an INI file. The INI file shown in Listing 1 represents a simulation that is performed on the Example_Network topology (network) with 100 qubits for each QNIC (buffers), and all qubits have some amount of measurement errors (Measurement_error_rate). The success probability of photon emission from one quantum memory is 0.8 (emission_success_probability). The paths are generated randomly (TrafficPattern), and end nodes measure 8000 Bell pairs for tomography at the end of Rule execution (distant_measure_count). This simulation stops after some amount of time in the simulation (sim-time-limit). Other than these variables, the default values (for example no errors on the link) are used for the simulation.
Once a user compiles the source, the simulation starts. All nodes are booted with the above parameters. QRSA is started, and if indicated tomography on the links (used to populate link information for routing) begins. For application traffic, each end node acts as an Initiator, randomly selects one other end node to be the Responder, and issues a "Connection Setup Request". For example, for the path highlighted in red in Figure 8, the Initiator issues a Connection Setup Request and sends it to the Responder.
The simulation automatically stops and invokes postprocessing functions to make statistics. An example output file is shown in Listing 2. Note that the fidelity in the listing is obtained from a density matrix constructed via quantum state tomography run by the application, and therefore corresponds to what real-world software would actually determine from its tests.

Correctness and performance
Based on the previous discussions on the simulator's design and usage, we now investigate its validity and performance. First, we give a simple numerical experiment to check whether the result out of the simulator agrees with theoretical values. After that, we extend the scale of the experiment to measure the performance.

Correctness
In order to test the correctness of QuISP we set up a small linear network composed of two end nodes A and C, with a single repeater B as depicted in Figure 1. The end nodes are separated from the repeater by equal distance d.
To make this test analytically tractable, we consider a simplified error model for the quantum link by assuming only Pauli X errors may affect the flying qubits with probability P X per kilometer of the fiber. This error model may transform the ideal Bell pair |Φ + = (|00 + |11 )/ √ 2 into a new state |Ψ + = (|01 + |10 )/ √ 2 and vice versa. The state of the two qubits after the flying qubit reaches repeater B is given by a statistical mixture of |Φ + with probability p clean , and |Ψ + with probability 1 − p clean , where Upon arrival of the flying qubits to the repeater node B they are measured together in the Bell basis. This results in the desired state |Φ + between qubits A-C if both pairs A-B 1 and B 2 -C are clean or if both suffered from an error [58]. Therefore the fidelity of qubit pair A-C is given by We simulated this network with P x = 0.02 and 7000 measurements. We repeated the simulation ten times. The resulting mean empirical fidelity and its standard deviation are shown in Figure 9 alongside the theoretical fidelity F AC . We varied the distance d from 1 to 10 kilometers. The empirical result matches the theoretical prediction.

Scalability
In order to show that QuISP can scale up to the goals mentioned in Section 4.1, we investigated the performance of QuISP in terms of events processed per second and the duration of CPU time taken to generate one end-to-end Bell pair. These results can indicate the relation between the scale of experiments and the time it takes to run the simulation. Although the total number of events depends on multiple factors such as the complexity of the protocol, network size, and number of buffer qubits each node has, it can give us some indications regarding the total work required.
For these tests, we decided to use the Docker environment QuISP provides. Table 2 details the simulation environment.
The work required to run the simulation can be broken down into three main parts: routine background task for link generation, connection setup, and RuleSet execution. The total simulation duration will be the aggregation of time required for each connection to establish, and each connection is scaled by number of hops, number of Bell pairs required, and the success probability of link generation. Since the load of RuleSet ... execution is directly related to the number of connections for any network size and connection setup is a very lightweight task, we decided to show the performance of QuISP using the simplest case, a linear topology with one connection. Figure 10 is the network that we use for performance analysis. We performed quantum state tomography using 200 E2E Bell pairs between these two end nodes with the probability of link generation at 0.32. We can adjust the number of repeaters between two end nodes, with the E2E distance growing as we add links. Each repeater has two QNICs with n Q qubits per QNIC. The end nodes contain one QNIC with n Q qubits. The total number of qubits is therefore N = 2n Q (n R + 1), where n R is the number of repeaters.
The RuleSet structure we used in this benchmark is similar to the one depicted in Figure 4. The Bell pairs are purified once before undergoing entanglement swapping at every step. The order of the entanglement swapping was chosen to be the binary tree structure [8]. The purification is also done on the end-to-end Bell pair shared between the two end nodes before tomography. Figure 11 shows the average CPU time (in seconds) per end-to-end Bell pair generated. When the number of repeaters increases, the amount of time to generate one end-to-end Bell pair increases because it includes far more operations (more entanglement swapping operations and more rounds of purification) than that with the smaller number of repeaters. In theory, the scaling of work required for generating end-to-end Bell pair with N repeaters is linearithmic (O(N log N)). We can clearly see that the scaling of our simulation time grows no worse than polynomially.
We also investigated the number of events processed per second. Figure 12 shows the number of events per second for different network sizes and number of qubits inside QNICs. Increasing the number of repeaters results in longer simulation time in the scaling as we expected. On the other hand, doubling the buffer qubits inside a QNIC (from 100 to 200) resulted in around twice the simulation time even though the total number of processed events are similar. This indicates that QuISP might have some kind of unintended overhead which scales linearly on the number of buffer qubits, which we expect to fix in a near-term release.

Conclusion
This paper proposed a discrete event simulator for a largescale quantum network on top of OMNeT++. QuISP is based on two essential aspects, the RuleSet protocol and error basis simulation. The RuleSet protocol allows us to establish end-to-end connections that operate fully asynchronously and autonomously. Error basis is a new representation of the quantum state, which only tracks the transition of errors rather than the full quantum state, enabling simulation that scales efficiently in both qubits involved in a single quantum state and in number of network nodes, links and connections.
We provided evidence that suggests QuISP properly models the real world by comparing simulation results with analytic results for a small network. We also assessed the performance of QuISP in terms of the actual CPU time and the number of events handled per second. While this cost varies with configuration, simulation cost does not row explosively, and simulations finish fast enough for researchers to explore aspects of the design space smoothly.
The simulator is a work in progress, but is complete enough for use as a research tool, and indeed our group has several projects running in parallel, each of which begins by defining the Rules that govern a new type of network functionality. Our current advanced projects are extending the functionality to include 2G, all-optical, multi-party protocols and internetworking.
Planned improvements can be split into three categories: core functionality, performance, and extensions. Connection teardown is incomplete, and routing, while correctly handling heterogeneity, is static; these limit our ability to compare multiplexing schemes effectively, and will be remedied soon. Performance measurements suggest substantial room to improve, and we expect to reach our goal of simulating an internet of 100 networks of 100 nodes each.
All of this is, of course, in pursuit of a functional, robust, extensible, and above all useful Quantum Internet. To get there, we must design and test architecture and protocols with emphasis on scale and heterogeneity. We expect that QuISP will contribute to advancing this design and implementation.

Acknowledgement
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA2386-19-1-4038.

Availability
The simulator source code, documentation and sample configurations are fully available online as an open-source project. Feedback, requests for features, and code contributions are welcome. Our mission is not only to build an advanced simulator, but also to make it easy to use for quantum networking research and as a teaching and learning tool. We are also alpha testing a web browser-based simulation interface created by compiling OMNeT++ into Web Assembly (wasm), allowing new users to learn about quantum networking without dealing with the complex build environment of OMNeT++.  Listing 3 is an example of NED syntax that defines the network with several quantum devices. Here, we have 11 EndNode, 5 Router, 6 Repeater and 8 HoM initialized with vectors of each module. The connections between these components are defined under connections which specify the port to be used and the type of channel such as ClassicalChannel to connect them.

A NED (Network Definition) in OMNeT++
The devices involved in the network are also defined as NED style language by nesting them as submodules and passing the parameters they require.In this examaple, address, nodeType and distance are given by network. The atomic components such as BSA, StationaryQubit are defined as simple, and the combination of simple modules such as QNode, Router are defined as module.

B Transition matrix Q
The transition matrix Q allows us to evolve the error probability vector in time and is given by Q =           P I P X P Y P Z P R P E P L P X P I P Z P Y P R P E P L P Y P Z P I P X P R P E P L P Z P Y P X P I P R P E P L 0 0 0 0 P I + P Z + P R P X + P Y + P E P L 0 0 0 0 P X + P Y + P R P I + P Z + P E P L 0 0 where P j represents the error probability per unit time. It captures the nature in which quantum errors interact with each other. For example, a qubit affected by a Pauli X error can be transformed into a clean qubit by another Pauli X error acting on it. This is reminiscent of how two classical bit flip errors cancel each other. Other errors may interact in more complicated ways. For example, a qubit that has been affected by the energy relaxation error R will stay in the same state if no further error occurs, or a Pauli Z occurs, or a the relaxation process takes place again. This is captured by the transition matrix element Q 55 = P I + P Z + P R . Some error states cannot be reached from certain error states. In particular, a qubit that has decohered via the relaxation or excitation process cannot go back to being affected by a simple Pauli error. This is captured by Q i j = 0. In this sense, the most destructive type of error is photon loss L captured by the fact that Q 7 j = 0 except for Q 77 = 1. Once a qubit is lost it we cannot recover it.


extfieldsofstudy:
Medicine
Computer Science

s2fieldsofstudy:
Computer Science
Medicine

text:
Weakly-supervised High-resolution Segmentation of Mammography Images for Breast Cancer Diagnosis

In the last few years, deep learning classifiers have shown promising results in image-based medical diagnosis. However, interpreting the outputs of these models remains a challenge. In cancer diagnosis, interpretability can be achieved by localizing the region of the input image responsible for the output, i.e. the location of a lesion. Alternatively, segmentation or detection models can be trained with pixel-wise annotations indicating the locations of malignant lesions. Unfortunately, acquiring such labels is labor-intensive and requires medical expertise. To overcome this difficulty, weakly-supervised localization can be utilized. These methods allow neural network classifiers to output saliency maps highlighting the regions of the input most relevant to the classification task (e.g. malignant lesions in mammograms) using only image-level labels (e.g. whether the patient has cancer or not) during training. When applied to high-resolution images, existing methods produce low-resolution saliency maps. This is problematic in applications in which suspicious lesions are small in relation to the image size. In this work, we introduce a novel neural network architecture to perform weakly-supervised segmentation of high-resolution images. The proposed model selects regions of interest via coarse-level localization, and then performs fine-grained segmentation of those regions. We apply this model to breast cancer diagnosis with screening mammography, and validate it on a large clinically-realistic dataset. Measured by Dice similarity score, our approach outperforms existing methods by a large margin in terms of localization performance of benign and malignant lesions, relatively improving the performance by 39.6% and 20.0%, respectively. Code and the weights of some of the models are available at https://github.com/nyukat/GLAM

Introduction
Convolutional neural networks (CNNs) have revolutionized medical image analysis (Bekkers et al., 2018;Ching et al., 2018;Topol, 2019;Prevedello et al., 2019;Lutnick et al., 2019;Karimi et al., 2020;Gaonkar et al., 2021). These networks achieve impressive results, but their outputs are often difficult to interpret, which is problematic for clinical decision making (Yao et al., 2018). Designing Figure 1: Comparison of saliency maps generated by CAM (Zhou et al., 2016), GMIC (Shen et al., 2021), and the proposed method on a mammography image containing a malignant lesion (first row, red) and a benign lesion (second row, green). Both CAM and GMIC produce coarse saliency maps that fail to localize the lesion accurately. The proposed method generates a high-resolution saliency map that precisely localizes the lesions.
In this work, we propose GLAM (Global-Local Activation Maps), a novel framework to generate fine-grained segmentation using only image-level labels. The proposed model processes high resolution medical images in a memory-efficient way. The main idea behind GLAM is to select informative regions (patches) that may contain ROIs via coarse-level localization and then to perform segmentation on selected patches rather than the entire image in a weakly supervised manner. We train and evaluate GLAM on a dataset containing more than one million mammography images. We demonstrate that the model outperforms existing baselines in segmentation of both benign and malignant lesions, improving the Dice similarity score relatively by 39.6% and 20%, respectively, while preserving classification accuracy. To achieve that, GLAM produces fine-grained saliency maps with a 300 times higher resolution than previous works (Shen et al., 2021) (736 × 480 pixels for 2944 × 1920 pixels input images). In Figure 1, we illustrate how the saliency maps generated by GLAM enable high-resolution segmentation of lesions relevant to breast cancer diagnosis.

Background
WSL is the task of learning to locate ROIs (i.e. the objects) in an image when only image-level labels are available during training. WSL methods are usually based on CNNs that produce saliency maps encoding the location of ROIs. To train the whole system using only image-level labels, the saliency maps are collapsed to predictions indicating the presence of each class using a pooling function. Once the CNN is trained, the saliency map can be used for localization (Oquab et al., 2015;Zhou et al., 2016). WSL has been applied in a wide range of medical-imaging applications, including the detection of lung disease in chest X-ray images (Wang et al., 2017;Yao et al., 2018;Tang et al., 2018;Ma et al., 2019;Liu et al., 2019;Guan et al., 2018), diagnosis of injuries from pelvic X-ray images (Wang et al., 2019), brain lesion segmentation (Wu et al., 2019a), breast MRI analysis (Luo et al., 2019), cancer detection in lung CT images (Feng et al., 2017;Schlemper et al., 2018), and scan-plane detection in ultrasound (Schlemper et al., 2018;Baumgartner et al., 2017).
The majority of these works focus on images that have relatively low resolution, that is, 512 × 512 pixels or less. Only a few works have considered higher resolution images, which are standard in some imaging procedures such as screening mammography (Shen et al., 2019(Shen et al., , 2021. In this work we focus on the diagnosis of breast cancer from screening mammography images. Breast cancer is the second leading cause of cancer-related deaths among women (Bray et al., 2018) and screening mammography is the main tool for its early detection (Marmot et al., 2013). CNN classifiers have shown promise in diagnosis from mammograms (Zhu et al., 2017;Kim et al., 2018;Ribli et al., 2018;Wu et al., 2019b;Geras et al., 2019;McKinney et al., 2020;Shen et al., 2019Shen et al., , 2021. Accurate localization of suspicious lesions is crucial to aid clinicians in interpreting model outputs, and can provide guidance for future diagnostic procedures. However, existing methods that explain their predictions, e.g. Shen et al. (2019Shen et al. ( , 2021, offer only coarse localization. GLAM is inspired by recent works (Yao et al., 2018;Shen et al., 2019Shen et al., , 2021Shamout et al., 2020), which improve classification accuracy by processing image patches selected from coarse saliency maps. The main innovation of GLAM with respect to these works is that it generates a a high-resolution saliency map from the selected patches, which significantly improves lesion segmentation accuracy.

Proposed Approach
Our goal is to generate fine-grained saliency maps that localize objects of interest in high-resolution images using only image-level labels during training. We start this section by describing the inference pipeline of our approach in Section 3.1. We then describe each module in detail in Section 3.2. Finally, we explain the training strategy in Section 3.3.

Inference pipeline
As illustrated in Figure 2, during inference, our system processes an input x ∈ R H,W as follows: 1. The image x is fed into the global module, a memory-efficient CNN denoted by f g , to produce an image-level coarse saliency map S g and an image-level class predictionŷ g . Figure 2: Inference pipeline of GLAM. 1) The global network f g is applied to the whole image x to obtain a coarse image-level segmentation map S g . 2) Based on this coarse-level segmentation, several patches are extracted from the input image.
3) The local network f l processes these patches to generate a high-resolution saliency map S l .
2. We select M patches from x based on S g . To do that, we greedily select the patches for which the sum of the entries in S g is the largest (see Algorithm 1 for a detailed description).
3. We feed the selected patchesx 1 , . . . ,x M to the local module f l , another CNN which produces a fine-grained saliency map associated with each patch. We then remap the patch-level saliency maps back to their location in the original input image. We denote the saliency map obtained through this procedure by S l . GLAM produces an image-level saliency map S g and a fine-grained multi-patch saliency map S l . These maps are aggregated through averaging to produce the final saliency map S c = (S g + S l )/2. In addition, GLAM generates a classification output, which is produced by the global module (as this yields the best classification accuracy).

Module parameterization
Global module The architecture of f g is based on the design of Shen et al. (2019Shen et al. ( , 2021, which is similar to ResNet (He et al., 2016) with a reduced number of channels for memory efficiency. The main difference between f g and the global module of Shen et al. (2019Shen et al. ( , 2021 is that we combine saliency maps at different scales to generate the global saliency map, inspired by Sedai et al. (2018) who observed that using convolutional feature maps extracted only from the last layer of a CNN may be suboptimal in localization of small objects. We generate saliency maps at different scales (S 0 , S 1 and S 2 ) using a pyramidal hierarchy of feature maps (see Figure 6 in Appendix B). Our image-level saliency map S g is obtained by averaging them. For each S n (n ∈ {0, 1, 2}), we obtain a classification predictionỹ n associated with S n using top t% pooling (see Appendix B), where t is a hyperparameter. The image-level classification predictionŷ g is calculated by averagingỹ 0 ,ỹ 1 , andỹ 2 . Additionally, we output a representation vector z g to feed it into a fusion module (described below) to enable joint training with the local module. See Appendix B for more details.
Local module Our local module is based on ResNet-34 with a reduced stride in the residual blocks to maintain a higher resolution. We replace the global average pooling and the last fully connected layer by a 1 × 1 convolution layer followed by a sigmoid non-linearity. The local module is applied to each of the selected patchesx k (k ∈ 1, . . . , K) to extract a patch-level saliency map A k . As we only have image-level labels, we need to train the patch-level network without patchlevel labels. To address this challenge, we combine insights from weakly-supervised localization and multiple-instance learning to train the patch-level saliency maps hierarchically. In multiple  instance learning (Maron and Lozano-Pérez, 1998) the labels are associated with bags of instances (the label is negative if all instances in a bag are negative, and positive otherwise). In our case each instance is a patch, and the patches from an image form a bag. We use a patch aggregation function f p to combine the information across all patches and form an image-level predictionŷ l . Formally, we haveŷ l = f p (A 1 , . . . , A k ). We propose two different patch aggregation functions.
• Attention-based aggregation: top t% pooling is applied to A k to obtain a patch-level prediction y k . Additionally, we output a representation vector z k for each patch, which will be aggregated to z l and fed into the fusion module. We use the Gated Attention Mechanism ((Ilse et al., 2018)) to combine the prediction and the representation vectors using attention weights α i ∈ [0, 1]. The prediction is given byŷ l = K i=1 α iŷi , and the representation vector by z l = K i=1 α i z i . Section 3.3 explains when we use each aggregation method. Refer to Appendix C for more details.
Fusion module In order to jointly optimize the parameters of the global and local modules, we incorporate a fusion module consisting of one fully connected layer that combines the representation vectors from the global (z g ) and local modules (z l ) to produce a fusion predictionŷ f . Formally, , where w f is a vector of learnable parameters.

Training strategy
The training strategy that achieves the best experimental performance for GLAM is sequential. We first train the global module, then we train the local module and finally we train both together. This makes sense because in order to train the local module f l effectively, we need to select meaningful input patches. Since the selection relies on the global saliency map, this requires pretraining the global module. Our training strategy is as follows (see also Figure 3).
1. Train the global module with the loss function L g = n∈{0,1,2} (BCE(y,ỹ n ) + λ (i,j) |S n (i, j)|) using the whole training set. Here, BCE(y,ỹ n ) is the binary cross-entropy loss and (i,j) |S n (i, j)| is an 1 regularization enforcing sparsity of the saliency map. The hyperparameter λ balances the two terms. The model is selected based on segmentation performance on the validation set.
2. Create training data for the local module. To produce negative examples, we select K patches randomly from images with negative labels. The positive examples are generated by selecting K patches from images with positive labels using Algorithm 1 applied to the saliency maps produced by the pretrained global module.
3. Train the local module minimizing the cost function L l = BCE(y,ŷ l ) + λ (i,j) |S l (i, j)|. Here we apply the concatenation-based aggregation described in Section 3.2 because it yields the best results experimentally and speeds up training.
4. Train the global and local modules jointly with the fusion module. Here we use the attentionbased aggregation in the local module. The fusion module outputs a fusion predictionŷ f . The loss function for joint training is Note that the number of patches selected as input to the local module during training (K) and inference (M ) do not need to be the same. In fact, we find empirically that it is beneficial to use K as large as possible (within the constraints imposed by GPU memory), which is consistent with Shen et al. (2021). However, choosing a large M increases the false positive rate. We used K = 6 and M = 1, based on experiments reported in Appendix G.4 and G.5.

Dataset and evaluation metrics
We trained and evaluated the proposed model on the NYU Breast Cancer Screening Dataset v1.0 (Wu et al., 2019c) that includes 229,426 exams (1,001,093 images) from 141,472 patients. Each exam contains at least four images with a resolution of 2944×1920 pixels, corresponding to the four standard views used in screening mammography: R-CC (right craniocaudal), L-CC (left craniocaudal), R-MLO (right mediolateral oblique) and L-MLO (left mediolateral oblique). The dataset is divided into disjoint training (186,816), validation (28,462) and test (14,148) sets, ensuring that each patient only belongs to one of the sets. Each breast has two binary labels indicating whether malignant or benign lesions are present. A subset of the images with lesions have pixel-level annotations provided by radiologists, which indicate the position of the lesions. Note that the dataset (458,852 breasts) contains more exams without lesions (452,311 compared to 5,556 with benign lesions, and 985 with malignant lesions). To account for this imbalance when training our models, at each epoch we use all exams with lesions and an equal number of randomly-sampled exams without lesions.
To measure classification performance, we report the area under the ROC curve (AUC) for identifying breasts with both malignant and benign lesions. To evaluate localization ability, we use the Dice similarity coefficient and pixel average precision (PxAP) (Choe et al., 2020). PxAP is the average of the area under the precision-recall curve for each pixel. The threshold to compute the precision and recall is either chosen for each image (image-level PxAP) or fixed for the whole test set (dataset-level PxAP). These metrics are described in more detail in Appendix E.

Comparison to baselines
We compare GLAM to two baselines: GMIC (Shen et al., 2021), a WSOL method specifically designed for high-resolution medical images, and CAM (Zhou et al., 2016), one of the most popular WSOL methods for natural images. We use the same backbone architecture for the two baselines and the global module of GLAM. The same hyperparameter tuning is applied to all models to ensure a fair comparison, as described in Appendix F. The models are selected based on segmentation performance, which is not necessarily equivalent to selection based on classification performance (see Appendix D for further discussion). We also compare to a U-Net (Ronneberger et al., 2015) trained with strong supervision using the pixel-level annotations. This provides an upper bound for the segmentation performance of the WSL methods. In Table 1, we report the performance of GLAM and the baselines. GLAM outperforms both GMIC and CAM in all segmentation evaluation metrics by a large margin, while achieving very similar classification accuracy. We observe a performance gap between GLAM and the strongly supervised U-Net model trained with ground-truth segmentation annotations, which indicates that there is room for further improvement.

Ablation study
We analyze different design choices in GLAM through an extensive ablation study. Due to space constraints, here we only discuss the advantages of training the global and local modules jointly, and the segmentation properties of the global and local saliency maps. We defer results on the following choices to the appendices: design of the global module (Appendix G.1), selection of training data for the local module (Appendix G.2), design of the local module (Appendix G.3), number of input patches used by the local module during training (Appendix G.4) and inference (Appendix G.5), the fusion module (Appendix G.6), and hyper-parameter selection: (Appendix G.7, G.8 and G.9).
Joint training of local and global modules. Here we compare the GLAM saliency map S c−joint obtained via joint training as described in Section 3.3 with (1) the global saliency map S g from the global module pretrained in isolation, (2) the local saliency map S l for the local module trained using patches selected using a frozen global module, (3) a saliency map S c−sep obtained by averaging S g and S l . The results reported in Table 2 show that S c−sep is superior to S g and S l , but joint training achieves better performance across all metrics.
Segmentation properties of the global and local saliency maps In Figure 4, we plot the Dice scores of the global (S g ) and local (S l ) saliency maps generated by GLAM for 400 randomly selected examples from the validation set. Each saliency map has different strengths and weaknesses. S l fails completely for a subset of examples, this is because the patch-selection procedure did not select the correct patches (lower row in Figure 5). On the remaining examples, S l tends to outperform S g because it has a much higher resolution. However, in some cases it underperforms S g , often in cases where the ground-truth segmentation is larger than the size of the patch (upper row in Figure 5). Averaging S l and S g (our strategy of choice in GLAM) achieves high-resolution segmentation, while hedging against the failure of the local module.

Conclusion
In this work, we propose a novel framework to perform weakly-supervised segmentation of images with very high resolution, which outperforms existing methods. Our results suggest that the general principle underlying GLAM (hierarchical selection of saliency maps) could be effective in other applications involving high-resolution images and videos. Another interesting question for future research is how to extend this principle to settings where high-resolution segmentation is desired, but (in contrast to breast cancer) the regions of interest do not tend to be small. 3. Additional details on the local module in Appendix C.
4. A discussion on the relationship between the localization and classification performance in Appendix D.

A detailed description of the evaluation metrics in Appendix E
6. A description of the hyperparameter tuning procedure in Appendix F.

A description of additional ablation studies in Appendix G.
Appendix A. Patch selection algorithm Algorithm 1 Patch selection algorithm: x denotes the input image, S denotes the saliency map that the patch selection method is based on (in our case of interest, this is the saliency map S g from the global module).x k denotes the selected patches cropped from the original image. K is the predefined number of selected patches.
Require: We perform patch selection using the same greedy algorithm (Algorithm 1) as in GMIC (Shen et al., 2021). In each iteration, we select the rectangular region in the coarse saliency map generated by the global module, which has the largest average intensity (see line 7 of Algorithm 1). Then the region is interpolated so that it maps to the corresponding location on the input image. Line 12 ensures that the extracted ROI patches do not significantly overlap with each other. It is worth noting that we do not assume that all patches extracted from a positive example contain tumors: only some of the K patches extracted by the global module will. This is accounted for by the patch aggregation function in the local module, which produces a single classification output for the K aggregated patches. A problem arises if the patch-selection algorithm fails to select any patches containing lesions for a positive example. This is why it is beneficial to use a larger number of patches during training (see Appendix H.4). It is worth noting, however, that this situation does not arise often in our experiments. In order to verify this, we used a set of images with lesions from the validation set, and checked the three first patches extracted by the global module. These patches did not contain any lesions only in 40 out of 300 images with benign lesions (13%) and in 6 out of 75 images with malignant lesions (8%).

Appendix B. Global module
Top t% pooling  Figure 6: Architecture of the global module. Both high-level (h 2 ∈ R 46×30×256 ) and low-level feature maps (h 0 ∈ R 184×120×64 , h 1 ∈ R 92×60×128 ) are utilized to obtain multi-scale saliency maps (S 0 ∈ R 184×120×2 , S 1 ∈ R 92×60×2 and S 2 ∈ R 46×30×2 ). For each scale n ∈ {0, 1, 2}, we use top t% pooling to transform saliency maps S n into class predictionsỹ n . A combined cross-entropy loss from the three scales is used for training. We compute the global saliency map S g by combining the individual saliency maps from three scales (S 0 , S 1 and S 2 ). We apply maximum pooling on the spatial dimension of h 2 ∈ R 46×30×256 to obtain the representation vector z g ∈ R 256 , which is fed to the fusion module during joint training.
As shown in Figure 6, the global module provides a multi-scale pyramidal hierarchy of feature maps (h 0 , h 1 , h 2 ) when processing an input image. The feature map corresponding to the deepest layer (h 2 ) has low resolution due to downsampling and therefore has limited localization ability. However, it is useful for classification, because it aggregates information from the whole image. Sedai et al. (2018) found that using only the deepest feature maps can negatively influence the localization of small objects. This motivates combining feature maps at several layers (in our case h 0 , h 1 , h 2 ) to produce saliency maps (Feng et al., 2017;Sedai et al., 2018;Yao et al., 2018). To do this, we use a training loss that includes classification estimates obtained from the three feature maps (see Eq. (1) below).
As shown in Figure 6, given a gray-scale input image of x ∈ R 2944×1920 pixels, we have feature maps in three scales (h 0 ∈ R 184×120×64 , h 1 ∈ R 92×60×128 and h 2 ∈ R 46×30×256 ). For each feature map h n of scale n, we use 1 × 1 convolution followed by a sigmoid function to transform h n to a saliency map S n ∈ [0, 1] h×w×|C| , where C = {malignant, benign}. Each element of S n , S n c (i, j), represents the contribution of a spatial location (i, j) towards classifying the input as class c ∈ C. We apply top t% pooling (Shen et al., 2021) as an aggregation function f agg (S n c ) : R h,w → [0, 1] to transform S n c to the image-level class prediction for class c. Formally, y c n = f agg (S n c ) = 1 where H + denotes the set containing the locations of the top t% values in S n c , where t is a hyperparameter. During training we minimize a sum of the cross-entropy losses between the image-level label y andỹ n for the three scales (n ∈ {0, 1, 2}). The training loss function of the global module is the following: (1) Here, BCE(y,ỹ n ) is the cross-entropy loss, (i,j) |S n (i, j)| is a 1 regularization term enforcing sparsity of the saliency map (Shen et al., 2021) and λ is a hyperparameter to balance the two loss terms.
We obtain the global module predictionŷ g by taking the average of the class prediction across the levelsŷ g = (ỹ 0 +ỹ 1 +ỹ 2 )/3. We obtain the global module saliency map S g by combining S 0 , S 1 and S 2 . Note that the resolutions of S 0 , S 1 and S 2 are different (S 0 ∈ R 184×120×2 , S 1 ∈ R 92×60×2 and S 2 ∈ R 46×30×2 ). In order to combine them, we upsample S 0 and S 1 to match the resolution of the S 2 using nearest-neighbor interpolation. Then, we compute the global module saliency map as S g = γ 0 S 0 + γ 1 S 1 + γ 2 S 2 . Here γ 0 , γ 1 and γ 2 are hyperparameters that can be tuned on the validation set. They should satisfy the conditions: γ 0 +γ 1 +γ 2 = 1 and γ 0 ≥ 0, γ 1 ≥ 0, γ 2 ≥ 0. We set γ 0 = 0.2, γ 1 = 0.6, γ 3 = 0.2 according to the Dice score on the validation set. We apply maximum pooling on the spatial dimension of h 2 ∈ R 46×30×256 to obtain a representation vector for the global module z g ∈ R 256 , which is fed into the fusion module during joint training of the global and local modules.
Utilizing multi-scale feature maps improves the segmentation performance of S g . As a result, we obtain a better feature map S c . This is demonstrated empirically in the ablation study reported in Appendix G.1.

Appendix C. Local module
We design the local module in order to preserve as much spatial resolution as possible. The resolution of the global saliency map is only 184 × 120 pixels. Therefore, the global module backbone is not suitable as a local module backbone. Typical classification networks also do not suit our needs. In our case, the size of the patches that are processed by the local module is 512 × 512. For such input resolution, the size of the saliency map produced by ResNet-34 is only 28 × 28 pixels. In order to preserve spatial resolution, we build the backbone for the local module by reducing the stride in all ResNet blocks of the ResNet-34 network to 1. When operating on an input patch with size 512 × 512 pixels, the resulting network outputs a saliency map with a resolution of 127 × 127 pixels. Note that we use ResNet-34 rather than ResNet-18 because it has a larger receptive field, which compensates for the reduced receptive field due to the reduced stride. We provide an ablation study of different choices for the local-module backbone architecture in Appendix G.3.
Architecture of the backbone of the local module is described in the left panel of Figure 7. In the right panel, we show two patch-aggregation strategies. We report an ablation study on the value  Figure 7: Architecture of the local module. The backbone network (left) is applied to each of the selected patchesx k (k ∈ 1, . . . , K) to extract a patch-level saliency map A k and a feature map h k . The patch-level saliency maps can be combined using two aggregation strategies: (1) Concatenation-based aggregation (top right), where we concatenate the saliency maps spatially and apply top t% pooling.
(2) Attention-based aggregation (bottom right), where top t% pooling is used to obtain patch-level predictionsŷ k from each patch, which are then combined using attention weights α i ∈ [0, 1] computed by the gated attention mechanism (Ilse et al., 2018). That is, the classification prediction is computed asŷ l = K i=1 α iŷi , and the representation vector as of t in the top t% pooling step in Appendix G.7. The selected value for both aggregation strategies is 20.

Appendix D. Localization and classification
A few works (Choe et al., 2020;Choe and Shim, 2019;Singh and Lee, 2017;Yao et al., 2018) observed that localization and classification accuracy of weakly-supervised models are not perfectly correlated. An empirical study by Choe et al. (2020) suggested that the best localization performance is usually achieved at early epochs of training, when the classifier is not fully trained. Furthermore, it has been observed that classification performance does not necessarily correlate with the localization performance across different architectures (Choe et al., 2020). We made similar observations in our work. Our local module reduces the stride of the convolution kernel to preserve more spatial resolution, which results in a smaller receptive field and has some adversarial effects for classification performance. Therefore, only the global module classification prediction is used as the classification result for our system. We provide the segmentation and classification results of different local module architectures in Appendix G.3. In this ablation study, we investigate the influence of the design of the local module on its segmentation and classification performance. We compare Resnet-34, a ResNet-18 with reduced stride in the residual blocks (ResNet-18-HR) and the proposed local module, which consists of a Resnet-34 with reduced stride in the residual blocks (Resnet-34-HR). Recall that the motivation to reduce the stride is to produce a saliency map with higher resolution, as explained in Appendix C. The results are shown in the following table. We see that carefully selecting the backbone architecture is crucial to ensure a good segmentation performance of the local module. It is worth mentioning that reducing the stride decreases the classification performance of the local module (which is why we do not use its output to perform classification). In this ablation study we use the model trained with six patches and evaluate its performance when using different numbers of patches during inference. The scatter plot in Figure 8 compares the segmentation performance when using one patch and when using three patches for 400 example images from the validation set. Using one patch leads to better performance for most of the images, but completely fails for some. Failure may occur when the patch selected by the global module for the positive example does not contain lesions. For our dataset, using one patch produces better segmentation performance overall, but this may not be the case for other applications (e.g. if there are a large number of lesions in each image).

G.6. Fusion module
In this ablation study, we compare the proposed model to two different models without fusion modules. In Model 1 the local module is trained using concatenation-based aggregation, in Model 2 it is trained with attention-based aggregation (see Section 3.2). To ensure a fair comparison, the models are obtained by jointly training the same pretrained global and local modules. We report the results in the following  Dice score (benign) Figure 8: Dice score of the local module when using one or three patches as its input during inference. The results correspond to 400 examples from the validation set. Using one patch leads to better performance for most of the images, but completely fails for a some. Failure may occur when the patch selected by the global module for the positive example does not contain lesions.

G.7. Impact of top t% pooling
As already shown by Shen et al. (2021), the choice of pooling hyperparameter has a significant impact on localization performance. In order to study how this hyperparameter influences the localization ability of the local module, we experimented with five choices of top t% pooling: t ∈ {10, 15, 20, 30, 60}. For all of them, we use the same trained global module, and train the local module without joint training. In the following table, we report the segmentation score achieved by each value of t on the test set. The results indicate that the pooling hyperparameter has a significant impact on the localization performance of the patch-level network (i.e. the local module), and the optimal t may be different for different classes. To investigate the influence of the λ hyperparameter associated with 1 regularization, we conduct the following experiments for the global module. We fixed all other hyper-parameter to the same values as used in the optimal model. To keep the comparison scenario simple, the saliency map of the global module is computed as S g = S 0 +S 1 +S 2

3
. λ is randomly selected from 10 [−5.5,−3.5] on a logarithmic scale during hyper-parameter searching for training, we take the following samples of λ for ablation study by randomly selection within the searching space: 10 −3.00 , 10 −3.53 , 10 −3.81 , 10 −4.24 , 10 −4.73 and report their Dice score on the test set. From the following table, we see that the performance of the saliency map is stable within the range of our hyper-parameter search space. For simplicity, in GLAM we set S c = (S g + S l )/2. Here, we investigate applying a weighted average S c = γ c S g + (1 − γ c )S l , with an additional hyperparameter γ c ∈ [0, 1]. The results in the following table show that this may slightly improve performance for some metrics.