Gustavo Girão

Daniel Barcelos
email: danielb@inf.ufrgs.br

Rech Flávio
email: flavio@inf.ufrgs.br

Wagner

Performance and Energy Evaluation of Memory Organizations in NoC-based MPSoCs under Latency and Task Migration

Keywords: Multiprocessor-System-on-Chip, Network-on-chip, Memory Organization, Cache Coherence, Task Migration, Performance and Energy Evaluation

This chapter presents a study on the performance and energy consumption arising from distinct memory organizations in an NoC-based MPSoC environment. This evaluation considers three sets of experiments. The first one evaluates the performance and energy efficiency of four different memory organizations in a situation where a single application is executed. In the second experiment, a traffic generator is responsible for the injection of synthetic traffic into the system, simulating the impact of the parallel execution of additional applications and increasing the latency of the NoC. Results show that, with a low NoC latency, the distributed memory presents better results for applications with low amount of data to be transferred. On the other hand, results suggest that shared and distributed shared memories present the best results for applications with high data transferring needs. In the second set of experiments, with higher NoC latency, for applications with low communication bandwidth requirements, a memory organization that is physically centralized and logically shared (called nDMA) is shown to have a smooth performance degradation when additional traffic rises up to 20% of the network capacity (22% degradation for an application demanding high communication, and 34% degradation for a low communication one). In contrast, a distributed memory model presents 2% of degradation in an application with high communication requirements, when traffic rises up to 20% of the network capacity, and reaches 19% of degradation in low communication ones. Shared and distributed shared memory models are shown to present lower tolerance to high latencies. A third set of experiments evaluates the performance of the four memory organization models in a situation of task migration, when a new application is launched and its tasks must be distributed among several nodes. Results show that the shared memory and distributed shared memory models have a better performance and energy savings than the distributed memory model in this situation. In addition, the nDMA memory model presents a smaller overhead when compared to the shared memory models and tends to reduce the traffic in the migration process due to the concentration of all memory modules in a single node of the network.

Introduction

Nowadays, embedded systems have become very complex. This complexity has many reasons, but the most evident one is the use of such devices for general purpose computing, leading to the execution of many different and complex applications. However, even with higher performance requirements, low power design is still a very desirable goal in portable devices [START_REF] Marwedel | Embedded System Design[END_REF].

To support processing requirements and also meet stringent constraints in terms of area and memory, as well as low energy consumption and low power dissipation, a solution using several cores in a single chip is widely adopted. This architecture is known as Multiprocessor System-on-Chip (MPSoC). This scenario usually implies a communication bandwidth between cores that demands a more efficient communication mechanism than a single bus [START_REF] Lee | On-Chip Communication Architechture Exploration: a Quantitative Exploration of Point-to-Point, Bus and Networkon-chip Architectures[END_REF]. With this concern in mind, the concept of Network-on-Chip (NoC) has been created.

Considering an MPSoC scenario, memory organization plays a key role since it is not only a major performance bottleneck but also represents a significant component in terms of energy consumption. In addition, memory organization is closely related to the communication model adopted in the application development. For instance, when using a shared memory organization, the communication mechanism usually adopted is the memory itself and, therefore, the memory organization becomes even more important.

Realizing that NoCs are communication structures with high scalability, it is not hard to imagine a situation with dozens or hundreds of processing elements and memory nodes, running a large number of applications concurrently. In this scenario, it is of great interest the evaluation of the behavior of different memory organizations when the network latency increases due to the large number of components and applications in the system. In addition, the memory model also impacts the system performance when a new application is dynamically launched and a task migration mechanism is applied such that a new task allocation is found which better meets system requirements, especially real-time and energy constraints.

This chapter presents a study on the performance and energy consumption arising from distinct memory organizations in an NoC-based MPSoC environment. This evaluation considers three sets of experiments, running on a virtual platform. The first one evaluates the performance and energy efficiency of four different memory organizations in a situation where a single application is executed. In the second set of experiments, a traffic generator is responsible for the injection of synthetic traffic into the system, simulating the impact of the parallel execution of additional applications and increasing the latency of the NoC.

The following memory organizations have been implemented in the virtual platform and evaluated in the experiments: (i) distributed memory, where processors have their local private memories; (ii) shared memory, with a single memory component in a dedicated node on the NoC that is accessed by all processors; (iii) distributed shared memory, composed by several physically distributed memory nodes that share the same address space; and, finally, (iv) a physically shared but logically distributed memory, whose communication model resembles a DMA communication protocol and is thus called nDMA.

Experiments show that, considering the communication requirements of an application, the results of performance and energy consumption may widely vary. For applications with high communication demands, the distributed memory model presents the highest tolerance to communication latency in most situations. However, for applications with low communication workload, the distributed memory model seems to present a larger degradation when NoC latency increases. Also, the nDMA model presents better results as the communication workload decreases. Experiments with the distributed memory model present a variation from 2% to 19% of performance reduction when traffic load rises from 10% to 20%, when using applications with high communication workload and low communication workload, respectively. On the other hand, the nDMA model shows a variation of performance reduction from 22% to 33% in the same situations, which represents a lower relative degradation if compared to the distributed memory model. On the other hand, the shared and distributed shared memory models present low tolerance to high latencies due to the use of a remote memory for communication among tasks.

A third set of experiments show that the shared memory and distributed shared memory models have a better performance and energy savings than the distributed memory model in a task migration situation. In addition, the nDMA memory model presents a smaller overhead when compared to the shared memory models and tends to reduce the traffic in the migration process due to the concentration of all memory modules in a single node of the network.

The remaining of this chapter is organized as follows. Section 2 discusses related work. Section 3 presents the virtual platform used to implement the experiments. Section 4 presents the experimental setup. In Sections 5 and 6, results for experiments with and without additional synthetic traffic, respectively, are presented. Finally, Section 7 draws conclusions and introduces future work.

Related Work

Several works regarding memory hierarchy in multiprocessor systems have been developed. However, the majority of these works only consider the use of busses instead of NoCs as communication mechanisms. In such systems, the massive communication parallelism may lead to different side effects due to the memory hierarchy.

Marescaux et al. [START_REF] Marescaux | The Impact of Higher Communication Layers on NoC Supported MPSoCs[END_REF] show a comparison between caches and scratchpads in an NoC-based MPSoC using a distributed shared memory model. In this environment, the use of six DSPs with local L1 caches and a shared L2 cache is considered. Experiments with two NoCs with different QoS methods are presented. The results show that scratchpads have a better performance than caches. However, these experiments do not consider a cache coherence mechanism in hardware. In this case, the adoption of a software coherence conservative approach that invalidates shared data on every access might have led to unnecessary invalidations.

Monchinero et al. [START_REF] Monchiero | Exploration of Distributed Shared Memory Architectures for NoC-based Multiprocessors[END_REF] explore the use of a distributed shared memory in an NoCbased MPSoC. The platform presents private L1 caches for each core and a shared L2 cache. Each processor has its own address space but there are also several banks of a distributed shared memory. A hardware MMU manages the shared data allocated to each processor. It is shown that by increasing the number of distributed shared banks the performance also increases but only up to some point where the NoC size leads to a greater latency. It is also shown that the energy consumption drops as the number of distributed shared banks increases.

Enright-Jerger et al. [START_REF] Enright-Jerger | Virtual Tree Coherence: Leveraging Regions and In-Network Multicast Trees for Scalable Cache Coherence[END_REF] propose a new cache coherence solution for multi-core architectures. The Virtual Tree Coherence (VTC) relies on a virtually ordered interconnection that keeps track of sharers of a coarse grain region. The protocol works in such a way that a virtual tree of nodes that share some region is established and each access to this shared region by any of these nodes leads to a message request to the root node of this virtual tree. The root node requests the data to the node that currently owns them. This request is performed as a multicast message, similarly to conventional snoop implementations. However, those multicast message requests are performed in a tree-like fashion in order to decrease the latency on the network, when compared to a multicast based on sequential unicast messages. The VTC solution is compared to a directory-based protocol and to a greedy-order protocol extended onto an unordered interconnect. VTC presented results 25% and 11% better, respectively, concerning performance. Nonetheless, this work does not present results about energy consumption of those cache coherence solutions.

Although the works presented in this section deal with memory organizations in an NoC-based MPSoC environment, none of them considers a high latency situation.

The SIMPLE Virtual Platform

Aiming at an accurate evaluation of the tolerance of memory models to a high latency scenario, four distinct memory models were implemented in the SIMPLE (Simple Multiprocessor Platform Environment) virtual platform. SIMPLE is a SystemC, cycle-accurate virtual platform that emulates an NoC-based MPSoC.

In SIMPLE, each Processing Element (PE) is a multi-cycle Java processor that is a hardware implementation of the Java Virtual Machine [START_REF] Ito | Making Java Work for Microcontroller Applications[END_REF]. Each instruction takes from 3 to 14 cycles (not including a possible cache miss). To generate the Java bytecodes, a compiler that follows the JVM specification is used. This compiler generates the contents of both the instruction memory and data memory customized for the application. These memories are used as inputs for the simulation in SIMPLE.

The NoC used in SIMPLE [START_REF] Zeferino | RASoC: a Router Soft-core for Networks-on-Chip[END_REF] implements a wormhole packet switching to reduce energy consumption. It also uses XY routing to avoid deadlock situations and a handshake control flow. Additionally, each router has five bi-directional ports with input buffer size of four phits. The phit size is four bytes.

As already mentioned, the simulator supports distinct memory organizations (distributed memory, shared memory, distributed shared memory, and nDMA memory) and cache configurations, regarding size, replacement policy, associativity, and block size.

For the distributed memory model, each router of the NoC is attached to a PE that, in turn, is attached to private memories (instruction and data). In this configuration, there are no caches and the communication mechanism uses message passing. Figure 1a depicts an example of this configuration.

The second model is a shared memory. Here, each router is attached to a PE or to a global Data Memory (the placement of each resource in the network is also configurable). Each PE has its own private data cache, while instruction memories are still local. In such environment there is a cache coherence problem. To solve this problem, SIMPLE adopts a directory-based cache coherence solution [START_REF] Girão | Cache Coherency Communication Cost in a NoC-based MPSoC Platform[END_REF]. This solution centralizes the memory access requests on an entity (the directory) that, based on the state of the block (if it is dirty or clean), makes decisions that could lead to invalidation or write-back requests if the block is dirty. The flowcharts for read and write operations are presented in Figures 2 and3, respectively. The shared memory in SIMPLE has a hardware-implemented directory, and each request from a cache to the memory is responded by it. This shared memory environment is shown in Figure 1b. The third memory configuration is a distributed shared memory, represented in Figure 1c. In this situation, there can be more than one shared data memory module. However, all of these modules share the same address space. This means that, if a memory module ends with the address N, then some other module in the system begins with the address N + 1. In this scenario, there are also private data caches for each PE and all global memory modules have their own directory module to maintain coherence. For each router on the NoC there is a Network Interface. In routers with a PE there is also a component called Memory Access Handler (MAH) that receives the signal from the processor when it needs data. The MAH module acts as an interface between the cache and the network interface when there is a cache miss. It creates a message to be sent to the directory module requesting a block. When the data are ready they are sent to the PE by the MAH module.

In both shared and distributed shared memory organizations the communication between the PEs is performed by means of shared variables protected by mutexes. The operations down and up in a mutex are made through test-and-set and test-andreset operations, respectively, supported by the hardware (caches and directory).

The fourth memory model available in SIMPLE is the nDMA organization, which implements a physically shared and logically distributed memory model, as depicted in Figure 1d. In this model the data memory is placed in an exclusive node, as in the shared memory model. This memory node, however, has N banks, each of them storing the memory data of each of the N PEs. Each memory bank has its own address space, and no PE can access the data of another PE. In addition, each PE has private data caches. When there is a miss, the cache sends a request to the memory node, where a Memory Controller (MC) receives the message, identifies the PEs that is requesting data, and, based on that, accesses the memory bank that hold its memory data.

A direct consequence of having different address spaces is that the most intuitive communication method to be used is message passing. Message passing is essentially implemented sending data from a processor (from its own memory, to be more precise) to another one through the physical communication mechanism. However, in the nDMA model, all memories are centralized in the same node, and, therefore, there is no need to send data through the network. Hence, the message passing method must be modified to work in such environment.

Basically, the processor that wishes to send a message must send a copy request message to the memory (instead of sending it to the target processor) informing an address source and the amount of bytes to be copied. In the memory node, the memory controller identifies this message and sends a message to the target processor informing the intention of the source processor to send data to it. The target processor replies with the address destination for that specific communication. When this reply message arrives at the memory node, the memory controller starts copying from one memory bank (containing the data memory of the source processor) to another one (containing the data memory of the target processor). Figure 4 illustrates this message passing mechanism. Since the memory node needs to be, at some level, programmed by the source processor to copy bytes from one address to another, the memory controller works similarly to a DMA, and, therefore, this memory organization is called nDMA (standing for NoC DMA).

Considering that each processor has its own private data cache, at the moment of sending the copy request the source processor must perform a forced write-back operation of the blocks inside the range of data that it is trying to send. However, the cache performs this write-back operation only on the blocks that have been modified. In a similar way, the target processor must also perform forced write-back operations on the blocks inside the range of the target addresses for that message, before the data are exchanged. Furthermore, the cache of the target processor must invalidate the blocks modified by the data exchange. These are the only moments during the entire communication process when data are exchanged through the network.

Experimental setup

The experiments consider four applications: a matrix multiplication, a motion estimation algorithm, a Mergesort algorithm, and a JPEG encoder.

The Matrix Multiplication was parallelized in such a way that each processor multiplies a subset of lines of matrix A by a subset of lines of matrix B. Each matrix used in simulations has 32 x 32 elements.

In the Motion Estimation, every PE searches a macroblock (a subset of an image) in a different part of the reference image. In the simulations, a macroblock of 8x8 pixels and an image in QCIF format (176x144 pixels) have been used.

For the Mergesort, the parallelism took advantage of its divide-and-conquer nature. Initially, each PE performs the Mergesort on a subset of the vector. Afterwards, one PE is responsible for assembling the whole vector, using the subsets already ordered by the various PEs.

A JPEG encoder can be seen as a three step algorithm. The first two steps (2-D DCT and Quantization) can be performed in parallel for different parts of the image. However, the third step (Entropy coding) can only be correctly performed with the whole image. Based on that, this parallel approach divides an image of 32x16 pixels in eight 8x8 blocks, and each PE is responsible for executing the first two steps on an equal amount of those eight blocks. At the end of those steps, each PE sends the resulting blocks to a master PE, which performs the final step with the complete image.

Considering the data inputs for the applications described above, Figure 5 represents their communication workload regarding different numbers of processors. Based on this chart, it is expected that the Motion Estimation algorithm will generate a larger amount of data exchanges.

Three kinds of experiments were performed. The first one evaluates the different memory organizations on an environment executing only a single parallelized application, with no addition of synthetic traffic. The second set of experiments investigates the behavior of the same memory organizations with the addition of different synthetic traffics, thus simulating the execution of several concurrent applications. Experiments have been performed for different numbers of processors and for different cache sizes. The third set evaluates the overhead in the system caused by a task migration. The task migration mechanism adopted in the experiments is quite simple. It is performed in three distinct steps that are performed sequentially by the processor were the task is located. The three steps correspond to the transmission (and proper receiving) of program code and data memory contents. Of course, the transmission of the data memory is not required when the shared data model is used. Also, because of the fact that the experiments assume a scenario where the tasks are recently created, there is no meaningful stack contents and therefore no need to migrate it. Experiments evaluate two characteristics: performance, measured by the total execution time of each application, and the overall dynamic energy spent, including processors, network, and memories. For the energy of the processors, a cycle-accurate power simulator [START_REF] Beck Filho | CACO-PS: a General purpose Cycle-accurate Configurable Power Simulator[END_REF] is used. For the network (including buffers, arbiter, crossbar, and links), the Orion library [START_REF] Wang | Orion: a Power-Performance Simulator for Interconnection Networks[END_REF] is applied, and for the memory and caches the Cacti tool [START_REF] Wilton | Cacti: An Enhanced Cache Access and Cycle Time Model[END_REF] is used. In all simulations, the processors and the NoC operate at 100 Mhz and the technology considered was 0.18µm. At this technology node, the static energy consumption is negligible. For applications with low communication workload (such as Matrix Multiplication, Mergesort, and JPEG), the distributed memory model presents better results than all other models. As for the Motion Estimation algorithm, the chart shown in Figure 6 indicates that the distributed memory does not present an overall result better than the other organizations. This is due to the fact that this algorithm requires a large number of data exchanges between processors, as depicted in Figure 5.

Experiments without Synthetic Traffic

The nDMA memory presents very similar results if compared to the other two shared memory models. For the Mergesort and JPEG application (Figures 8 and9, respectively), the nDMA presents even better results, especially considering more realistic MPSoC sizes, as with eight PEs. Considering the other applications (like the Matrix Multiplication depicted in Figure 7), the distributed shared memory shows a slightly superior performance, especially with eight processors. This enhanced performance is due to the parallel memory accesses performed by the several PEs. However, this advantage does not increase proportionally to the number of memory modules because not all of the modules have the same access profile. This means that some module may concentrate more accesses than others, and, hence, the performance speed up of this solution is not so high.

Energy Consumption

Concerning dynamic energy, Tables 1 thru 4 show the consumption for the four applications. For each application, there are four columns representing the energy consumption (considering an average between the cache sizes) for the Processors, Caches, Memories, and NoC for eight PEs.

Analyzing the energy consumption difference between the nDMA solution and the other two shared memory organizations, it is possible to see that the NoC consumption is much different. This behavior is a consequence of the fact that the nDMA memory uses small control messages for communication and does not need cache coherence messages. The use of a cache coherence solution in the shared and distributed shared memories also leads to invalidations that increase the memory accesses (in order to retrieve the block again) and the cache accesses. Therefore, the energy consumption of those components is also higher.

As depicted, in the distributed memory model the processor consumption is the major responsible for the energy consumption in the system. The NoC and sometimes also the memory energy are negligible in this environment. The memory energy is more significant in the Motion Estimation simulation due to the amount of data to be manipulated. For the two shared memory organizations and the nDMA model, the processor also plays a key role in the overall energy consumption. However, in the case of Matrix Multiplication, because of the larger caches, these components present a very significant energy consumption as well. In fact, for the shared memory organizations, the cache energy increases in a non-linear fashion as the cache size increases. As opposed to the increase of the cache sizes, the NoC energy decreases, and, except for the smallest cache size, one can say that the NoC is not a key factor in this case. The same can be said about the memory.

As an overall analysis, the energy consumption results tend to follow the performance results, and, therefore, the distributed memory presents better results for the low communication workload applications (Matrix Multiplication, Mergesort, and JPEG) and worst results for the Motion Estimation.

Experiments with Synthetic Traffic

The physical communication mechanism of a multiprocessor system directly affects the memory model adopted. This is due to the fact that the memory model leads to a certain communication model, which, in turn, is influenced by the physical communication mechanism. Based on that, this study tries to emulate an environment with high latency and quantify its impact on the system for different memory models.

In order to emulate such environment, a traffic generator was developed. The general idea is to create a synthetic traffic to increase the latency of the communication mechanism, thus emulating the parallel execution of several applications that generate communications among the various processors.

The Traffic Generator is placed inside the network interface (which is present in every node of the system), as depicted in Figure 10, and a state machine coordinates the inclusion of a traffic package in the Send Buffer. In addition, making use of an identifier in the header of the package, the traffic generator analyzes the incoming packages in the Receive Buffer to exclude any synthetic traffic package. This avoids the resource associated to the Network Interface (PE or a memory) from reading and processing it. The traffic generator woks in a very simple fashion. At every period (previously configured), it creates a 10-byte package to be sent to a destination node determined through a round-robin sequence of all routers in the system, one at a time. Each traffic package has a 10-byte size in order to guarantee that the network interface is able to send a package in only one cycle.

With this round-robin system, at every 10 cycles router 0 in the system sends a synthetic traffic packet to router 1, router 1 sends a packet to router 2, and so on. In the next period, router 0 sends a packet to router 2, router 1 sends a packet to router 3, and so on. Therefore, at a certain time, each router sends a synthetic traffic packet to a different destination. On the other hand, when these packets arrive at their final destinations, the Traffic Generator removes them from the Receive Buffer.

The experiments use 10% and 20% traffic loads, meaning that at every 10 cycles or 5 cycles, respectively, each router in the system sends a 10-byte traffic packet to a different router. Results indicate by how much the performance decreases (or the energy consumption increases), when the traffic load increases from 10% to 20%.

This Traffic Generator creates a traffic that is uniformly distributed both in time and in space. As future work, more complex generators (as in [START_REF] Carara | Communication Models in Networks-on-Chip[END_REF][START_REF] Mahadevan | A Network Traffic Generator Model for Fast Network-on-Chip Simulation[END_REF]) will be used, to assess the possible influence of the traffic model on the experimental results.

Performance

Considering the summarized results of reduction in performance depicted in Table 5, it is possible to see how the increase in the latency impacts the performance.

Although the distributed memory presents better results when compared to the other memory models, it is very clear how much the higher latency impacts the performance as the communication workload of an application increases. In the Motion Estimation application, a higher latency only affects the distributed memory results in 2%. However, when an application with less communication workload (the JPEG) is executed, the impact of a higher latency reaches up to 19%.

The shared memory and distributed shared memory models seem to suffer less from this communication workload variation, although, in absolute numbers, the overall decrease of performance is higher than in any other memory model. On the other hand, the nDMA model suffers less as the latency increases, and its absolute results are in the middle between the shared memory models and the distributed one. This is mainly due to the necessity of accessing a remote memory even though small control messages are used.

These results suggest that a higher latency in the NoC affects much more the distributed memory model than the other models, as the communication workload of an application decreases. Furthermore, the nDMA model seems to have results with less variation as the communication workload of an application decreases.

Energy Consumption

Again, the pattern present on the performance results also appears in the dynamic energy consumption. These results suggest that applications with high communication requirements seem to have a smaller reduction in the energy consumption as the latency increases.

In the same way as for the performance results, it is possible to see that a higher latency impacts the energy consumption of the memory models differently, depending on the communication workload characteristics of the various applications. Results in Table 6 show that the distributed memory model seems to suffer more in a high latency situation if the communication workload of the application is low. Again, the shared memory and distributed memory organizations seem to suffer more than other models, whereas the nDMA model presents a graceful degradation if the communication workload of the application is low. This section presents experimental results regarding the impact of the distinct memory models on task migration, considering the performance on the migration execution and the dynamic energy spent on it. Experiments quantify the task migration overhead in each model. Since all memory models consider that each PE has its own private instruction memory, all of them demand the transmission of program code during task migration, as well as of data memory, depending on the model. The experiments consider a worst case scenario, corresponding to the creation of a new application with N tasks in a single node and the following migration of these tasks to other nodes. Considering that N represents the number of processors, each processor has exactly one new task assigned to it after the migration. In this case, the new task has no data on the stack and therefore, it is not necessary to migrate its contents.

Performance results

According to the performance results presented in Figures 11 thru 14, the distributed memory model presents worst results in all cases. This was already expected due to the fact that the distributed memory model demands the migration of the whole data memory. This situation does not happen in the case of shared and distributed shared memories. In the case of the nDMA model, due to the fact that it is essentially a shared memory model, in a sense that all data (from all processors) is located in a single node, the copy can be made simultaneously to the program code migration. In addition, in the nDMA memory model there is no traffic overhead for the data memory transfer during migration. However, these are not the only factors responsible for the worse performance of the distributed memory model. Another factor is the size of the code, which, in the case of distributed memory, is usually higher than in other models. This is due to the fact that communication needs to be completely explicit in the code. The programmer has to describe the copying of data to sending buffers and the reading of data from receiving buffers, as well as calls for these functions. This procedure is less intense in the case of the nDMA memory model, which, although not requiring the transmission of the contents of the data memory, works by sending explicit messages to control communication between the processing elements. On the other hand, the shared memory and distributed shared memory models present the best performance in all cases. The fact that only the application code contents is copied is not the single reason for that, but also the simplicity of the code in the shared memory models.

Despite the better results of the shared and distributed shared memory models, it is important to note that the migration of the contents of the caches was not taken into account in these experiments. Therefore, a loss of performance may be expected during the initial stages of task execution after its migration in these cases, due to a considerable amount of compulsory cache misses. This case has not been considered in the experiments since the migration occurs at the time of the creation of the task and so the caches contain virtually no information relevant to that task. However, if one takes into account a migration during the execution of this task, the compulsory cache misses will occur. This is a situation whose outcome depends on the point in time within the execution of the task in which the decision is taken to perform the migration. This choice is a function of the task allocation algorithm and was not considered in the experiments, which evaluate only the mechanism of task migration itself.

Energy consumption results

Regarding the dynamic energy consumption, results shown in Figures 15 thru 18 demonstrate the same pattern presented in the previous performance results.

Among the components of the system, it is noticeable that the NoC has an extremely low energy consumption during the migration process when compared to the energy consumption of processors and memory. In particular, the migration process appears to be very expensive for the processor that needs to execute the whole process of sending and receiving messages and perform the routines of reading and writing on program and data memories (in the case of distributed memory).

Again, the results point to a worse result for the distributed memory, due to the reasons already discussed in the previous section. The shared and distributed shared memories are the most efficient ones, given the smaller size of code to be read and written in memory (thus reducing the number of memory accesses). This also reduces the execution time of the processor performing this task, which leads to a lower power consumption.

Conclusions and Future Work

This chapter presented a study on the tolerance of different memory organization models under high latency NoC scenarios and also the impact of these memory models in a task migration situation. The experiments considered not only different memory models and traffic loads but also different numbers of processors and cache sizes and various applications with different bandwidth requirements.

The experiments suggest that the distributed memory model presents a higher tolerance to communication latency in most situations. The nDMA model also presents some tolerance and presents better results as the communication workload increases. The high level of communication requirements on shared and distributed shared memory leads to the worst results as the synthetic traffic increases. This is due to the fact that not only network communication is required to retrieve application data but also the communication among tasks is performed through a remote memory.

Regarding dynamic energy consumption, it is possible to conclude that the most affected component of an MPSoC based on a high latency NoC is the NoC itself, while other components such as memory and cache are not affected, due to the fact that the number of accesses does not change and thus neither the dynamic energy consumption.

Although the distributed memory model presents better results in most situations, there is a pattern that shows that, for applications with low communication workload, it suffers more with higher latencies on the NoC. On the other hand, the nDMA memory model presents results with lower degradation as the communication workload of an application is low.

The experimental results also suggest that the distributed memory model is the one in which the task migration is more expensive. These experimental results were expected due to the simple fact that the distributed memory needs to send task data explicitly through the network, which is not true for the other models. In the case of shared and distributed shared memories the fact of having a global memory limits the amount of information to be transferred in a migration to the program memory. In the case of the model nDMA the situation is very similar. The difference occurs because the data memory copy is performed locally on the node where the memory banks are located, making it a quicker process, independent of the network. Another factor that explains why the shared and distributed shared memories obtain better results is the fact that the application code in such cases is simpler than in the distributed memory model, where all communication must be specified explicitly.

Future work includes more experiments considering other applications and synthetic traffics with more complex distributions, as well as the proposal and evaluation of an optimized task migration model for each memory organization, also considering high latency situations. The evaluation of the static energy consumption for more recent technology nodes must also be considered.

An ongoing work will propose a cluster scenario for an MPSoC with dozens or even hundreds of processors. According to the results presented in this chapter, it was possible to see that different applications have distinct memory hierarchy demands, and, therefore, an MPSoC with different memory organization solutions in the same chip, in order to address the needs of different applications running concurrently, makes sense. The ongoing work relies on these observations to build a system with multiple clusters on the NoC, where each cluster may have a distinct memory organization. Thus, it is possible to take advantage of parallelism at application and task level. Figure 19 gives an insight on how this MPSoC could be. In the figure, there are three applications allocated in the MPSoC and the processors running tasks from the same application build a single cluster. As the parallelism of an application changes during the execution [START_REF] Kumar | Single-ISA Heterogeneous Multi-Core Architectures: The Potential for Processor Power Reduction[END_REF], the task allocation can be dynamic, making the cluster size also dynamically modifiable, thus favoring another application that may need more processing or memory resources.

Each of these clusters makes use of a memory organization that provides a better performance for the task (or set of tasks) being performed on it. To accomplish the communication between clusters it is necessary to provide a component that builds a bridge between the memory models. This component could be centralized, so that all inter-cluster communication would be sent to a single component in the system. Another solution would provide components for inter-cluster communication within each cluster, in a manner similar to a router, but at a higher level of abstraction. These components would have to perform communication in accordance to the memory organization adopted by the cluster (message passing, shared variable, etc).

In addition, as future MPSoCs tend to provide adaptability, the memory hierarchy can be also dynamic in a sense that each memory node (i.e. a node that holds a general memory structure) can be configurable in terms of a distinct memory organization. Hence, the memory node can become a cache or a main memory, leading to a shared memory or to a distributed memory model, respectively. In order to do that, it is also required a software layer to identify the application needs.

Another objective of this future work is the support of consistency of data between different models of memory organization, as proposed in [START_REF] Dutt | Memory-aware NoC Exploration and Design[END_REF]. This problem can be seen as a more expanded version of the problem of consistent caches. However, in this case it is intended to maintain consistency between the memory hierarchies of several clusters. Again the solution to this problem is likely to be centralized, creating a template directory that keeps a record of the location of all data in the system. An alternative solution is the adoption of a hierarchical system of directories. In this case, each cluster contains a directory that centralizes data consistency for the cluster itself and responds directly to a central directory that maintains the consistency of data at a macro level.

Fig 1 .

 1 Fig 1. Memory models available in SIMPLE.

Fig. 2 .

 2 Fig. 2. Read operation using directory.

Fig. 3 .

 3 Fig. 3. Write operation using directory.

Fig. 4 .

 4 Fig. 4. nDMA message passing mechanism.

Fig. 5 .

 5 Fig. 5. Communication workload.

5. 1 Performance

 1 Concerning performance, Figures 6 thru 9 show the overall execution time of each application, measured in millions of cycles. For each memory organization, there are three columns representing different cache sizes -256, 512, and 1024 bytes in the Motion Estimation, Mergesort, and JPEG simulations, and 1024, 2048, and 4096 bytes for the Matrix Multiplication.

Fig. 6 .

 6 Fig. 6. Performance for Motion Estimation.

Fig. 7 .

 7 Fig. 7. Performance for Matrix Multiplication.

Fig. 8 .

 8 Fig. 8. Performance for Mergesort.

Fig. 9 .

 9 Fig. 9. Performance for JPEG.

Fig. 11 .

 11 Fig. 11. Task migration performance for the Motion Estimation.

Fig. 12 .

 12 Fig. 12. Task migration performance for the JPEG.

Fig. 13 .

 13 Fig. 13. Task migration performance for the Matrix Multiplication.

Fig. 14 .

 14 Fig. 14. Task migration performance for the Mergesort.

Fig. 15 .

 15 Fig. 15. Task migration energy consumption of the Motion Estimation.

Fig. 16 .

 16 Fig. 16. Task migration energy consumption of the JPEG.

Fig. 17 .

 17 Fig. 17. Task migration energy consumption of the Matrix Multiplication.

Fig. 18 .

 18 Fig. 18. Task migration energy consumption of the Mergesort.

Fig. 19 .

 19 Fig. 19. Cluster allocation in a large MPSoC.

Table 1 .

 1 Energy consumption (mJ) without synthetic traffic for Motion Estimation.

Table 2 .

 2 Energy consumption (mJ) without synthetic traffic for Motion Estimation.

Table 3 .

 3 Energy consumption (mJ) without synthetic traffic for Mergesort.

	Memory Model	µP	Motion Estimation D$ Mem NoC	Total
	Distributed Memory	3180 87%	--			452 12%	4 0.1%	3634 100%
	Shared Memory	668 79%	137 16%		13 2%	25 3%	843 100%
	Distr. Shared Memory	667 80%	136 16%		10 1%	24 3%	837 100%
	nDMA	866 75%	183 16%		96 8%	12 1%	1157 100%
	Memory Model	µP	Matrix Multiplication D$ Mem NoC	Total
	Distributed Memory	288 54%	--			243 45%	2 0.5%	533 100%
	Shared Memory	919 41%	1214 54%		21 1%	103 4%	2257 100%
	Distr. Shared Memory	919 41%	1214 54%	14 0.6%	83 4%	2230 100%
	nDMA	1342 44%	1628 54%		38 1%	24 1%	3032 100%
	Memory Model		µP	D$		Mergesort Mem NoC Total
	Distributed Memory		408 70%	--		173 29%	2 0.5%	583 100%
	Shared Memory		1260 54%	479 20%	122 5%	472 21%	2333 100%
	Distr. Shared Memory		1264 56%	484 21%	60 3%	440 20%	2244 100%
	nDMA		866 73%	238 20%	45 4%	42 3%	1191 100%

Table 4 .

 4 Energy consumption (mJ) without synthetic traffic for JPEG.

Table 5 .

 5 Impact of NoC latency on performance.

	Application	Motion Estimation Matrix Multiplication Mergesort JPEG
	Distributed Memory	2%	4%	6%	19%
	Shared Memory	29%	40%	40%	44%
	Dist. Shared Memory	33%	48%	49%	55%
	nDMA	22%	33%	33%	33%

Note: values in the table indicate how much the performance is reduced when traffic load increases from 10% to 20%.

Table 6 .

 6 Impact of NoC latency on energy consumption.

	Application	Motion Estimation Matrix Multiplication Mergesort JPEG
	Distributed Memory	1%	10%	14%	32%
	Shared Memory	25%	50%	65%	67%
	Dist. Shared Memory	44%	70%	77%	79%
	nDMA	15%	24%	24%	24%
	Note: values in the table indicate how much the energy consumption increases when traffic load increases
	from 10% to 20%.


Julian Krebs
email: julian.krebs@inria.fr

Tommaso Mansi

Hervé Delingette

Li Zhang

Florin C Ghesu

Shun Miao

Andreas Maier

Nicholas Ayache

Rui Liao

Ali Kamen

Robust non-rigid registration through agent-based action learning

Robust image registration in medical imaging is essential for comparison or fusion of images, acquired from various perspectives, modalities or at different times. Typically, an objective function needs to be minimized assuming specific a priori deformation models and predefined or learned similarity measures. However, these approaches have difficulties to cope with large deformations or a large variability in appearance. Using modern deep learning (DL) methods with automated feature design, these limitations could be resolved by learning the intrinsic mapping solely from experience. We investigate in this paper how DL could help organ-specific (ROI-specific) deformable registration, to solve motion compensation or atlas-based segmentation problems for instance in prostate diagnosis. An artificial agent is trained to solve the task of non-rigid registration by exploring the parametric space of a statistical deformation model built from training data. Since it is difficult to extract trustworthy ground-truth deformation fields, we present a training scheme with a large number of synthetically deformed image pairs requiring only a small number of real inter-subject pairs. Our approach was tested on inter-subject registration of prostate MR data and reached a median DICE score of .88 in 2-D and .76 in 3-D, therefore showing improved results compared to state-of-the-art registration algorithms.

Introduction

Registration of images with focus on the ROI is essential in fusion and atlasbased segmentation (e.g. [START_REF] Tian | A fully automatic multi-atlas based segmentation method for prostate MR images[END_REF]). Traditional algorithms try to compute the dense mapping between two images by minimizing an objective function with regard to some similarity criterion. However, besides challenges of solving the ill-posed and non-convex problem many approaches have difficulties in handling large deformations or large variability in appearance. Recently, promising results using deep representation learning have been presented for learning similarity metrics [START_REF] Simonovsky | A deep metric for multimodal registration[END_REF], predicting the optical flow [START_REF] Dosovitskiy | Flownet: Learning optical flow with convolutional networks[END_REF] or the large deformation diffeomorphic metric mapping-momentum [START_REF] Yang | Fast predictive image registration[END_REF]. These approaches either only partially remove the above-mentioned limitations as they stick to an energy minimization framework (cf. [START_REF] Simonovsky | A deep metric for multimodal registration[END_REF]) or rely on a large number of training samples derived from existing registration results (cf. [START_REF] Dosovitskiy | Flownet: Learning optical flow with convolutional networks[END_REF][START_REF] Yang | Fast predictive image registration[END_REF]).

Inspired by the recent works in reinforcement learning [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF][START_REF] Ghesu | An artificial agent for anatomical landmark detection in medical images[END_REF], we propose a reformulation of the non-rigid registration problem following a similar methodology as in 3-D rigid registration of [START_REF] Liao | An artificial agent for robust image registration[END_REF]: in order to optimize the parameters of a deformation model we apply an artificial agent -solely learned from experience -that does not require explicitly designed similarity measures, regularization and optimization strategy. Trained in a supervised way the agent explores the space of deformations by choosing from a set of actions that update the parameters. By iteratively selecting actions, the agent moves on a trajectory towards the final deformation parameters. To decide which action to take we present a deep dual-stream neural network for implicit image correspondence learning. This work generalizes [START_REF] Liao | An artificial agent for robust image registration[END_REF] to non-rigid registration problems by using a larger number of actions with a low-dimensional parametric deformation model. Since ground-truth (GT) deformation fields are typically not available for deformable registration and training based on landmark-aligned images as in rigid registration (cf. [START_REF] Liao | An artificial agent for robust image registration[END_REF]) is not applicable, we propose a novel GT generator combining synthetically deformed and real image pairs. The GT deformation parameters of the real training pairs were extracted by constraining existing registration algorithms with known correspondences in the ROI in order to get the best possible organ-focused results. Thus, the main contributions of this work are: [START_REF] Dosovitskiy | Flownet: Learning optical flow with convolutional networks[END_REF] The creation and use of a low-dimensional parametric statistical deformation model for organ-focused deep learning-based non-rigid registration. (2) A ground truth generator which allows generating millions of synthetically deformed training samples requiring only a few (<1000) real deformation estimations. (3) A novel way of fuzzy action control.

Method

Training Artificial Agents

Image registration consists in finding a spatial transformation T θ , parameterized by θ ∈ R d which best warps the moving image M as to match the fixed image F. Traditionally, this is done by minimizing an objective function of the form:

arg min θ F(θ, M, F) = D (F, M • T θ ) + R(T θ )
with the image similarity metric D and a regularizer R. In many cases, an iterative scheme is applied where at each iteration t the current parameter value θ t is updated through gradient descent: θ t+1 = θ t + λ∇F(θ t , M t , F) where M t is the deformed moving image at time step t: M • T θt .

Inspired by [START_REF] Liao | An artificial agent for robust image registration[END_REF], we propose an alternative approach to optimize θ based on an artificial agent which decides to perform a simple action a t at each iteration t consisting in applying a fixed increment δθ at : θ t+1 = θ t + δθ at . If θ is a ddimensional vector of parameters, we define 2d possible actions a ∈ A such that δθ 2i [j] = i δ j i and δθ 2i+1 [j] =i δ j i with i ∈ {0..d -1}. In other words the application of an action a t increases or decreases a specific parameter within θ t by a fixed amount where δ j i is an additional scaling factor per dimension that is set to 1 in our experiments but could be used e.g. to allow larger magnitudes first and smaller in later iterations for fine-tuning the registration.

The difficulty in this approach lies into selecting the action a t as function of the current state s t consisting of the fixed and current moving image: s t = (F, M t ). To this end, the framework models a Markov decision process (MDP), where the agent interacts with an environment getting feedbacks for each action. In reinforcement learning (RL) the best action is selected based on the maximization of the quality function a t = arg max a∈A Q (s t , a). In the most general setting, this optimal action-value function is computed based on the reward function defined between two states R(s 1 , a, s 2 ) which serves as the feed-back signal for the agent to quantify the improvement or worsening when applying a certain action. Thus, Q (s t , a) may take into account the immediate but also future rewards starting from state s t , as to evaluate the performance of an action a.

Recently, in RL powerful deep neural networks have been presented that approximate the optimal Q [6]. Ghesu et al. [START_REF] Ghesu | An artificial agent for anatomical landmark detection in medical images[END_REF] used deep reinforcement learning (DRL) for landmark detection in 2-D medical images. In the rigid registration approach by Liao et al. [START_REF] Liao | An artificial agent for robust image registration[END_REF] the agent's actions are defined as translation and rotation movements of the moving image in order to match the fixed image.

In this work, the quality function y a (s t ) ≈ Q (s t , a) is learned in a supervised manner through a deep regression network. More precisely, we adopt a single-stage MDP for which Q (s t , a) = R(s t , a, s t+1 ), implying that only the immediate reward, i.e. the next best action, is accounted for. During training, a batch of random states, pairs of F and M, is considered with known transformation T θ GT (with F ≈ M • T θ GT ). The target quality is defined such that actions that bring the parameters closer to its ground truth value are rewarded:

Q (s t , a) = R(s t , a, s t+1 ) = θ GT -θ st 2 -θ GT -θ a st+1 2 . (1)
The training loss function consists of the sum of L 2 -norms between the explicitly computed Q-values (Eq. 1) for all actions a ∈ A and the network's quality predictions y a (s t ) per action. Having a training batch B with random states s b the loss is defined as:

L = s b ∈B a∈A y a (s b ) -Q (s b , a)
2 . In testing, the agent iteratively selects the best action, updates the parameter θ t and warps the moving image M t as to converge to a final parameter set representing the best mapping from moving to fixed image (see Fig. 1b).

Statistical Deformation Model

One challenge of the proposed framework is to find a low dimensional representation of non-rigid transformations to minimize the number of possible actions (equal to 2d), while keeping enough degrees of freedom to correctly match images. In this work, we base our registration method on statistical deformation models (SDM) defined from Free Form Deformations (FFD). Other parametrizations could work as well. Typically, the dense displacement field is defined as the I1 𝒊 11

I2

...

Base Image

Intra-Subj.

Inter-Subj. summation of tensor products of cubic B -splines on a rectangular grid. Rueckert et al. [START_REF] Rueckert | Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration[END_REF] proposed to further reduce the dimensionality by constructing an SDM through a principal component analysis (PCA) on the B -spline displacements.

We propose to use the modes of the PCA as the parameter vector θ describing the transformation T θ that the agent aims to optimize. The agent's basic increment per action i is normalized according to the mean value of each mode estimated in training. To have a stochastic exploration of the parameter space, predicted actions a t are selected given fixed probabilities (see [START_REF] Liao | An artificial agent for robust image registration[END_REF]).

Fuzzy Action Control. Since parameters θ are the amplitudes of principal components, the deviation of θ 2m and θ 2m+1 from the mean µ m should stay within k-times the standard deviation σ m in testing. In order to keep θ inside this reasonable parametric space of the SDM, we propose fuzzy action controlling. Thus, actions that push parameter values of θ outside that space, are stochastically penalized -after being predicted by the network. Inspired by rejection sampling, if an action a moves parameter θ m to a value f m , then this move is accepted if a random number generated between [0, 1] is less than the ratio N (f m ; µ m , σ m )/N (h; µ m , σ m ) where h m = µ m + kσ m , and N is the Gaussian distribution function. Therefore, if |f m -µ m | ≤ kσ m , the ratio is greater than 1 and the action is accepted. If |f m -µ m | > kσ m then the action is randomly accepted, but with a decreased likelihood as f m moves far away from µ m . This stochastic thresholding is performed for all actions at each iteration and rejection is translated into adding a large negative value to the quality function y a . The factor k controls the tightness of the parametric space and is empirically chosen as 1.5. By introducing fuzzy action control, the MDP gets more robust since the agent's access to the less known subspace of the SDM is restricted.

Training Data Generation

Since it is difficult to get trustworthy ground-truth (GT) deformation parameters θ GT for training, we propose to generate two different kinds of training pairs: Inter-and intra-subject pairs where in both moving and fixed images are synthetically deformed. The latter pairs serve as a data augmentation method to improve the generalization of the neural network.

In order to produce the ground truth deformations of the available training images, one possibility would be to apply existing registration algorithms with optimally tuned parameters. However, this would imply that the trained artificial agent would only be as good as those already available algorithms. Instead, we make use of manually segmented regions of interest (ROI) available for both pairs of images. By constraining the registration algorithms to enforce the correspondence between the 2 ROIs (for instance by artificially outlining the ROIs in images as brighter voxels or using point correspondences in the ROI), the estimated registration improves significantly around the ROI. From the resulting deformations represented on an FFD grid, the d principal components are extracted. Finally, these modes are used to generate the synthetic training samples by warping the original training images based on randomly drawn deformation samples according to the SDM. Amplitudes of the modes are bounded to not exceed the variations experienced in the real image pairs, similar to [START_REF] Rueckert | Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration[END_REF].

Intra-subject training pairs can be all combinations of synthetically deformed images of the same subject. Since the ground-truth deformation parameters are exactly known, it is guaranteed that the agent learns correct deformations. In the case of inter-patient pairs a synthetic deformed image i mb of one subject I m is allowed to be paired with any synthetic deformed image i nc of any other subject I n with b, c denoting random synthetic deformations (see Fig. 1a). Thereby, the GT parameters θ GT for image pair (i mb , i nc ) are extracted via composition of the different known deformations such that ((i mb

• T i mb ,Im θ ) • T Im,In θ ) • T In,inc θ .
Note the first deformation would require the inverse of a known deformation that we approximate by its opposite parameters for reasons of computational efficiency. The additional error due to this approximation, computed on a few pairs, remained below 2% in terms of the DICE score.

Mini-batches are created online -during training -via random image pairing where intra-and inter-subject pairs are selected with the same probabilities. Through online random pairing the experience of new pairs is enforced since the number of possible image combinations can be extremely high (e.g. 10 12 ) depending on the number of synthetic deformations.

Experiments

We focused on organ-centered registration of MR prostate images in 2-D and 3-D with the use case of image fusion and atlas-based segmentation [START_REF] Tian | A fully automatic multi-atlas based segmentation method for prostate MR images[END_REF]. The task is very challenging since texture and anatomical appearance can vary a lot. 25 volumes were selected from the MICCAI challenge PROMISE124 and 16 from the Prostate-3T database5 including segmentations. Same images and the cases with rectal probes were excluded. Randomly 8 cases were chosen for testing (56 pairs), 33 for training. As preprocessing, translation-based registration for all pairs was carried out in 3-D using the elastix -framework [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF] with standard parameters followed by cropping and down sampling the images (to 100x100/75x75x20 pixels in 2-D/3-D respectively). For the 2-D experiments, the middle slice of each volume has been taken. For the purpose of GT generation mutual information as similarity metric and a bending energy metric was used. The optimization function was further constrained by a Euclidean point correspondence metric. Therefore, equally distributed points were extracted from the given mask surfaces. elastix was used to retrieve the solution with the weights 1, 3 and 0.2 for the above-mentioned metrics and a B-spline spacing of 16x16(x8) voxels. As a surrogate measure of registration performance we used the DICE score and Hausdorff distance (HD) on the prostate region. The extracted GT resulted in median DICE coefficients of .96 in 2-D and .88 in 3-D. Given the B-spline displacements, the PCA was trained with d = 15 modes in 2-D, d = 25 in 3-D (leading to 30 respectively 50 actions with a reconstruction error < 5% (DICE score) as a compromise to keep the number of modes relatively small. The network's two independent processing streams contained 3 convolutional (with 32, 64, 64 filters and kernel size 3) and 2 max-pooling layers for feature extraction. The concatenated outputs of the two streams were processed in 3 fully-connected layers (with 128, 128, 64 knots) resulting in an output with size 2d (equals the number of actions). Batch normalization and ReLu units were used in all layers. The mini-batch size was 65/30 (2-D/3-D). For updating the network weights, we used the adaptive learning rate gradient-based method RMSprop. The learning rate was 0.001 with a decay factor of 0.8 every 10k mini-batch back-propagations. Training took about 12 hours/ 1 day for 2-D and 3-D respectively. All experiments were implemented in Python using the deep learning library Theano including Lasagne6 . DL tasks ran on GPUs (NVIDIA GeForce GTX TITAN X ). During testing 200 MDP iterations (incl. resampling of the moving image) took 10 seconds (GPU) in 2-D and 90 seconds in 3-D (GPU). The number of testing steps was set empirically since registration results only change marginally when increasing the number of steps. In empirical 2-D experiments with 1000 steps the agent's convergence was observable.

For testing, the initial translation registration was done with elastix by registering each of the test images to an arbitrarily chosen template from the training base. Table 1 shows that our method reaches a median DICE coefficient of .88/.76 in 2-D/3-D and therefore shows similar performance as in [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF] with the best reported median DICE of .76 on a different data set. However, on our challenging test data our method outperformed the LCC-Demons [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF] algorithm with Table 1: Results of prostate MR registration on the 56 testing pairs. 2-D and 3-D results in comparison to elastix with B-spline spacing of 8 (e8) or 16 (e16) as proposed in [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF] and the LCC-Demons [START_REF] Lorenzi | LCC-Demons: a robust and accurate symmetric diffeomorphic registration algorithm[END_REF] algorithm (dem). T are the initial scores after translation registration with elastix. 3-D* are results with perfect rigid alignment T*. nfc are our results with no fuzzy action control (HD in mm). 4.3 6.8 6.7 2.5 3.9 10.7 10.9 3.9 2.9 2.3 6.8 6.4 2.5 2.2 1.9

manually tuned parameters and elastix, using similar parameters as proposed for prostate registration [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF] using B-spline spacing of 8 and 16 pixels. We found that better rigid registration can significantly improve the algorithm's performance as shown in the experiments with perfect rigid alignment according to the segmentation (3-D*). Extreme results are shown visually in Fig. 2.

Regarding the results of elastix and LCC-Demons, a rising DICE score was observed while HD increased due to local spikes introduced in the masks (visible in Fig. 2b) as we focused on the DICE scores during optimization for fair comparisons. In the 3-D* setting, DICE scores and HDs improved when applying fuzzy action control compared to not applying any constraints (see Table 1).

Conclusion

In this work, we presented a generic learning-based framework using an artificial agent for approaching organ-focused non-rigid registration tasks appearing in image fusion and atlas-based segmentation. The proposed method overcomes limitations of traditional algorithms by learning optimal features for decisionmaking. Therefore, segmentation or handcrafted features are not required for the registration during testing. Additionally, we propose a novel ground-truth generator to learn from synthetically deformed and inter-subject image pairs.

In conclusion, we evaluated our approach on inter-subject registration of prostate MR images showing first promising results in 2-D and 3-D. In future work, the deformation parametrization needs to be further evaluated. Rigid registration as in [START_REF] Liao | An artificial agent for robust image registration[END_REF] could be included in the network or applied as preprocessing to improve results as shown in the experiments. Besides, the extension to multimodal registration is desirable.

Disclaimer This feature is based on research and is not commercially available. Due to regulatory reasons its future availability cannot be guaranteed.

Fig. 1 :

 1 Fig. 1: (a) Training Data Generation: Synthetic (blue arrows) and intersubject GT deformations (black) are used for intra-(green) and inter-subject (red) image pairs for training. (b) Dual-stream network used for Q-value prediction ya including complete single-stage Markov Decision Process for testing (blue background).

(a) 2

 2 -D: Moving, Fixed, elastix -e8 (.84), elastix -e16 (.70), ours (.94). (b) 3-D: Moving, Fixed, elastix -e8 (.49), elastix -e16 (.59), LCC-Demons (.67), ours (.79).

Fig. 2 :

 2 Fig. 2: 2-D and 3-D registration results of extreme cases with segmentation masks overlays (fixed: green, moving: orange) and DICE scores in parenthesis.

  DICE Mean .76 .74 .77 .87 .62 .63 .64 .67 .75 .74 .72 .67 .79 .79 .80 DICE Med. .78 .79 .81 .88 .61 .71 70 .67 .76 .75 .77 .76 .80 .79 .81 DICE StD. .10 .15 .13 .05 .11 .22 .20 .11 .06 .08 .17 .23 .07 .05 .04 HD Mean 11.6 15.2 14.5 7.7 16.1 21.2 25.3 15.9 11.8 9.2 13.4 14.5 10.4 8.9 8.0

		2-D	3-D	3-D*
		T e16 e8 our T e16 e8 dem our T* e16 e8 dem nfc our
	HD Med.	11.7 13.2 13.0 7.2 15.2 18.0 21.7 15.8 11.2 9.0 11.6 13.5 10.8 8.8 7.9
	HD StD.

https://promise12.grand-challenge.org/

https://wiki.cancerimagingarchive.net/display/Public/Prostate-3T

https://lasagne.readthedocs.io/en/latest


Emmanuel Guillaud

Laurent Joly

Dominique De Ligny

Samy Merabia

Assessment of elastic models in supercooled water: A molecular dynamics study with the TIP4P/2005f force field

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Many materials display a viscoelastic behavior: they are elastic -i.e. they behave like solids -at short times, whereas they are viscous -i.e. they flow -at longer times. The distinction between liquid and solid is only conditioned by the timescale considered, the temperature and the level of stress applied to the material: under normal conditions, the Earth's mantle is solid at our timescale but flows on geological timescales, whereas water is liquid at our timescale and responds elastically at the picosecond scale [START_REF] Oswald | Rheophysics: the deformation and flow of matter[END_REF] . Glass former liquids are well known examples of such materials: at the laboratory timescale, they are glassy -i.e. solids -at low temperature, whereas they flow like viscous liquids above the glass transition temperature T g . Water falls in the category of glass formers, in the sense that it can be supercooled down to the homogeneous nucleation temperature T = 231 K 2 .

One of the most puzzling features of the glass transition is the huge slowing down of the dynamics as probed e.g. by the α-relaxation time or the shear viscosity. It has been proposed that this slowing down is strongly related to the short time elastic properties of the supercooled liquid [3][4][5] , the correlation being described by the so-called elastic models. In the case of supercooled water, its slow dynamics has been characterized either experimentally or through molecular simulations [6][7][8][9][START_REF] Kumar | Proc. Natl. Acad. Sci[END_REF][START_REF] Morón | [END_REF][12][13][14] . However, very few works addressed the connection between the slow water dynamics and short times elastic properties. Recently, Klameth and Vogel discussed a connection between the α-relaxation time and the high frequency shear modulus on the basis of Schweizer's elastic model 15 . We can also mention the study by Ciamarra et al., which does not concern specifically water, but a model glass former displaying a density maximum, as observed in water 16 .

The common idea behind elastic models is the fol-lowing [3][4][5] . Collective relaxation processes such as flow events are infrequent in the supercooled liquid and occur on long timescales, consisting in local molecular rearrangements characterized by an effective activation energy E a (T ). Elastic models relate the activation energy to the elastic properties of the supercooled liquid calculated on short -picosecond -time scales. In the shoving model for instance, the activation energy is proportional to the work W necessary to shove the surrounding of a given molecule, so that flow can occur. The work W depends on the elastic properties of the supercooled liquid, as described by the shear modulus, denoted G, and the bulk modulus, denoted K. If the local structure around a rearranging region has a spherical symmetry, the rearrangement consists in a purely shear deformation 17 . Consequently, the observables characterizing the slow dynamics, including the shear viscosity η or the collective structural relaxation time τ α -generally denoted X in the following, should be expressed as:

     X = X 0 exp E a (T ) k B T , E a (T ) = λ G p (T ) a 3 (1)
where k B denotes Boltzmann's constant, T is the temperature, a denotes a molecular length, G p (T ) is the high frequency plateau modulus and λ is a dimensionless number which takes values close to one (see Ref. 4 for a review of the various models leading to this expression). Another popular version of elastic models emphasizes correlations between the slow dynamics and the mean square displacement (Hall-Wolynes equation) 16,18 :

X = X 0 exp a 2 2 u 2 (T ) . ( 2
)
where u 2 (T ) is the vibrational mean square displacement. Elastic models, in their two versions, proved to provide a good description of various systems including metallic, polymer and colloid glasses 5,[19][20][21] , but showed also some limitations, as in the case of oxide glasses 22 .

In this paper, we assess the validity of the two versions of elastic models to describe the slow dynamics of supercooled water using molecular dynamics simulations. We consider the promising TIP4P/2005f force field which has been shown to reproduce accurately the water shear viscosity over a broad range of temperatures, extending deep in the supercooled regime 14 . We report a strong correlation between the observables characterizing the dynamics and the high frequency shear modulus or the mean square displacement for the temperatures considered in this study. Conversely, elastic models perform well at low temperatures, but for the viscosity two regimes with different effective molecular lengths a are observed, with a transition temperature on the order of 260 K. Interestingly, this temperature is close to the one where the Stokes-Einstein relation begins to be violated.

II. NUMERICAL METHODS

All calculations were performed using the LAMMPS package 23 .

Water molecules were modeled by the TIP4P/2005f force field 24 , which recently proved to be very precise to describe the shear viscosity and the selfdiffusion coefficient of supercooled water 14 . Cubic simulation boxes with periodic boundary conditions and containing respectively 8 3 , 10 3 , 12 3 , 14 3 , 16 3 , 18 3 , 20 3 and 22 3 molecules were equilibrated in the isothermal-isobaric ensemble under 1 atm pressure and various temperatures between 225 and 360 K, with four replica of the system for each box size. For the lowest temperatures considered, we did not detect any sign of ice nucleation. The selfdiffusion coefficient, the shear viscosity, the α-relaxation time and the high frequency shear modulus were then computed for each replica in the microcanonical ensemble. In detail, the self-diffusion coefficient was calculated from the slope of the molecular mean-squared displacement, the shear viscosity was computed within the Green-Kubo formalism, and the α-relaxation time was obtained from the coherent intermediate scattering function using a procedure detailed in Ref. 14. Table I summarizes the simulation times at different temperatures.

Finally, we extracted the high frequency shear modulus from the autocorrelation function of the traceless stress tensor components σ αβ (t) σ αβ (0) , where σ αβ denotes one of the five independent components of the traceless stress tensor (σ xy , σ yz , σ xz , (σ xx -σ yy )/2 and (σ yy -σ zz )/2, see Ref. 26). The plateau shear modulus G p is traditionally defined as the value of the elastic modulus where the derivative of log (G(t)) is minimum 5,27 . The corresponding plateau should be located around 170 fs, as inter-molecular connectivity bands appear around 200 cm -1 on the infrared absorption spectrum of bulk water at ambient temperature 28 . However, due to zero wave vector optical phonon modes 25 , the au- tocorrelation function of the stress tensor components strongly oscillates at short timescales (see Fig. 1), so that the plateau value G p must be extrapolated. We fitted the stress autocorrelation function for t > 0.03 ps

σ αβ (t) σ αβ (0) = G 0 k B T V exp - t τ M γ , (3)
where V denotes the volume of the system and T the temperature, which was found to give good results over long time scales for each temperature. In the following, we will therefore use the parameter G 0 as a proxy to the high frequency plateau modulus G p . Note that for the highest temperatures, extraction of G 0 becomes difficult as the oscillations of the stress tensor auto-correlation function extend over a time interval longer than the decay time τ M . The stretching exponent γ as a function of temperature is shown in Fig. 1. At high temperatures T > 280 K, the exponent γ is close to 1 and the stress autocorrelation function can be well described by a single exponential. Below 280 K, γ decreases with cooling, indicating that a single Maxwell relaxation time can not describe stress relaxation at low temperatures, as also observed in Refs 29,30. Details about the different methods used are provided in Ref. 14, where it is shown in particular that for viscosity measurements, the Green-Kubo formalism provides values in quantitative agreement with steady-state shear simulations.

III. RESULTS AND DISCUSSION

We first look at direct correlations between the shear viscosity or the α-relaxation time of liquid water with the high frequency shear modulus across Eq. ( 1). As we showed in a recent paper 14 , the TIP4P/2005f force field predicts a non-Arrhenius behavior for the shear viscosity and the α-relaxation time, even if deviations to the Arrhenius law are smaller in the case of α-relaxation time. Both observables are displayed in Fig. 2. The evolution of the high frequency shear modulus is also reported in the inset of this figure. Interestingly, the shear modulus G 0 decreases monotonously with increasing temperature. Therefore, the temperature dependence of G 0 follows that of the plateau modulus G p , and G 0 is clearly distinct from the infinite frequency shear modulus G ∞ 5,27 . Using the initial value G 0 of the stretched exponentials as a proxy of G p , the agreement between elastic models and our calculations is only partial: for the shear viscosity, Eq. ( 1) does not fit well the simulation data over the whole temperature range analyzed, and two temperature domains can be distinguished, as shown in Fig. 3. In the FIG. 2. Evolution of the α-relaxation time and the shear viscosity with respect to the temperature. Both observables display a non-Arrhenius behavior, although deviations to the Arrhenius law are smaller for the α-relaxation time. Inset: temperature dependence of the shear modulus G0.

FIG. 3. Correlations between the shear viscosity / αrelaxation time and the high frequency shear modulus of

water. While the α-relaxation time can be described with Eq. ( 1) over the whole temperature range, two different regimes are observed for the shear viscosity, with a transition temperature of ca. 261 K.

two domains, Eq. ( 1) describes the correlation between G 0 /T and the viscosity, but the parameter a is not the same at high and low temperatures. For T < 261 K, a ≈ 1.8 Å, while for T > 261 K, a ≈ 1.5 Å (assuming λ = 1). The fact that the parameter a is larger at low temperatures may be explained by an increase of cooperativity, as discussed in the following paragraph. As for α-relaxation time, Eq. ( 1) provides a good description of the slowing down of supercooled water with a single set of parameter a ≈ 1.6 Å over the entire temperature domain.

Interestingly, the temperature of transition between the high temperature and low temperature behavior of the viscosity is found to be close to 261 K, which corresponds to the temperature below which the self-diffusion coefficient D and the shear viscosity η decouple with respect to the Stokes-Einstein relation, Dη ∼ T , which we estimated to be 265 ± 15 K and related to the appearance of dynamical heterogeneities in a previous work 14 . The Stokes-Einstein relation describes a regime where particle motion is not cooperative and may be compared to a Brownian particle diffusing in a bath of neighbor molecules. Conversely, elastic models traditionally describe a situation where atomic motions involve collective displacements of the surrounding particles. However, here we saw that elastic models could still describe the slowing down of the dynamics at high temperatures, provided we defined an effective length a which takes a smaller value as compared with the low temperature value. The fact that a is smaller at high temperatures could then be related to the loss of cooperativity. Therefore, two regimes emerge in the description of the dynamics in the framework of elastic models: a diffusive homogeneous regime at high temperature where molecule displacements are almost independent and the Stokes-Einstein relation is obeyed, and a low temperature regime where molecule motions are cooperative.

We now assess the validity of the mean square displacement version of the elastic model, as described in Eq. 2. The quantity of interest here is the vibrational mean squared displacement, which can be accessed through the mean square displacement r 2 (t) , as illustrated in the top inset of Fig. 4. In a dense system, molecules move ballistically at short times while they diffuse at long times: r 2 (t) = 6Dt where D is the diffusion coefficient. At intermediate time scales, the mean square displacement of a molecule is limited by the presence of neighbouring molecules. As seen in fig. 4, at these intermediate time scales the MSD has a minimum slope, which corresponds to the time t DW when the particle shoves its neighbors. The corresponding value of the MSD is called mean square vibrational amplitude or Debye-Waller (DW) factor 5,17 , and denoted u 2 in the following.

Since the pioneering work of Buchenau and Zorn 31 , Eq. ( 2) has been shown to hold in many systems 16,32,33 . For supercooled water, Fig. 4 reports a good agreement between Eq. (2) and our simulation data at low temperatures. As for the previous elastic model, we find two different behaviors depending on the observable X. If X is the α-relaxation time, the agreement extends over the whole range of temperatures analyzed, with a molecular length a = 1.62 Å. In the case of the shear viscosity, two regimes with different molecular lengths appear (a = 1.97 Å at low temperature, and a = 1.48 Å FIG. 4. Correlations between the shear viscosity / αrelaxation time and the vibrational mean squared displacement. The vibrational mean squared displacement is determined as the value of the mean squared displacement of the center of mass of a water molecule when the slope is minimum (see top left inset). Similarly to the previous correlations with the elastic modulus, the α-relaxation time is described by Eq. ( 2) over the whole temperature range, while a transition is observed at ca. 258 K for the shear viscosity. The inset shows the correlation between the vibrational mean squared displacement and the high frequency shear modulus.

at high temperature), with a crossover temperature close to 258 K. This crossover temperature is very close to the one determined in the previous approach. The equivalence between the two approaches discussed is confirmed by the inset of Fig. 4, which demonstrates the good correlation between G 0 /T and the MSD. The difference of behavior between the shear viscosity and the α-relaxation time is not surprising, and can be related to the decoupling between these two observables observed at low temperatures 14,34 . The relevance of the vibrational MSD u 2 is supposed to hold while τ α is lower than t DW 32 , which corresponds to a temperature much higher than the crossover temperature T = 258 K. Actually, the activation energy associated with the two versions of the elastic models at this threshold temperature is close to the energy of a hydrogen bond (ca. 20 kJ/mol), which could explain the change of regime at high temperatures, since the connectivity of water molecules is mainly due to hydrogen bonds. At high temperature, thermal motion breaks the hydrogen bond network, so that the local tetragonal structure is short lived, and molecular displacements are no longer collective 30 .

Before concluding, we should say a word about the role of quantum effects in the dynamics of supercooled water. Indeed, our simulations are fully classical and one may ask the relevance of quantum effects to describe the slow dynamics of water. These include the description of the nuclear degrees of motion, the hydrogen bond net-work but also the temperature dependence of the mean square displacement. We have already shown that the TIP4P/2005f potential reproduces accurately the shear viscosity of supercooled water 14 . It would be interesting to measure the elastic properties of supercooled water at low temperatures, to see if the TIP4P/2005f potential can also describe the temperature dependence of the shear modulus, or the vibrational mean square displacement. This would allow us to conclude whether quantum effects do play a signicant role in the slow dynamics of water, or if they play a negligible role due to compensation effects. If quantum effects turn out to be important, demanding ab initio simulations could be unavoidable to investigate supercooled water dynamics [35][36][37] .

IV. CONCLUSION

In conclusion, we show here that the long time dynamics of liquid water at low temperatures is highly correlated with its short time elastic properties, quantified by the shear modulus at high frequency. In particular, it is surprising to see that in spite of the numerous anomalies of water, elastic models perform as well as for simple glass formers including metallic glasses or supercooled polymers. Moreover, for the viscosity two different temperature regimes with different effective molecular lengths a are observed, and the crossover temperature is close to the one where the Stokes-Einstein relation starts to be violated. These results call for experimental investigation of the high frequency shear modulus of supercooled water using e.g. inelastic X-ray scattering. Looking for correlations between the elastic modulus and the shear viscosity measured at low temperatures may also help in assessing the role of quantum effects which are not considered in our classical approach. Finally, measuring the high frequency elastic properties of water may shed light on the molecular mechanisms relevant to interfacial thermal transport at solid-water interfaces 38 .

FIG. 1 .

 1 FIG.1. Extraction of the high frequency shear modulus from the shear stress auto-correlation function. A stretched exponential fit describes data at long timescales. Top: exponent γ of the stretched exponential as a function of temperature. Middle: Shear stress auto-correlation function at high temperatures, T = 360 K and 300 K. The solid lines show the stretched exponential fits. Bottom: Shear stress autocorrelation function at T = 240 K, together with the stretched exponential fit. Inset: the stretched exponential functional form is compared to the MD data at short timescales, emphasizing the relatively large amplitude of the oscillations of the stress auto-correlation function due to zero wave vector optical phonon modes25 .

TABLE I .

 I Relaxation time τα and simulation times at different temperatures. Note that at each temperature, we considered eight simulation box sizes and four replicas of the system for each size.T (K) 225 240 255 270 285 300 315 330 345 360 τ α (ps) 31.8 18.3 6.27 3.04 2.08 1.08 0.748 0.498 0.413 0.376 simulation time (ps) 500 340 300 300 300 300 300 300 300 300 using a stretched exponential functional form,

ACKNOWLEDGMENTS

This study was performed using the computing resources of the PSMN, Computing Center of the Ecole Normale Supérieure in Lyon, and benefits from financial support of DFH/UFA 39 .