Eric Orellana-Romero

Javier Sanmartin-Hernandez

Cristian Duran-Faundez

Vincent Lecuire
email: vincent.lecuire@univ-lorraine.fr

Katherine Zapata-Quiñones

Evaluation of Block Interleaving Techniques for Robust Image Communication in Wireless Camera Sensor Networks

Keywords: Image communication, block interleaving, error resilience, wireless camera sensors, packet loss

We study block interleaving techniques for enhancing robustness against packet loss on wireless vision sensor networks. We evaluate three techniques: DSJ-AL interleaving, Turner and Peterson (TP) technique, and Torus Automorphisms (TA). Results in terms of image quality, obtained trough extensive simulations over monochromatic images, show the improvements provided by interleaving techniques in comparison with a sequential communication, with a slight favor to TA.

I. INTRODUCTION

Wireless Sensor Networks (WSNs) consist of a large number of small electronic devices with sensing, communication and processing capabilities [START_REF] Akyildiz | Wireless sensor networks: A survey[END_REF]. Usually, their application fields relate to simple scalar measures such as temperature, pressure or humidity [START_REF] Arampatzis | A survey of applications of wireless sensors and wireless sensor networks[END_REF], however, advances in microelectronics allow low cost capture and transmission of digital images, making possible Wireless Vision Sensor Networks (WVSNs) [START_REF] Soro | A survey of visual sensor networks[END_REF]. Such networks introduce a large number of open challenges, because the many problems present in traditional WSNs are increased by the use of more complex sensors (cameras) and the need to deal with the thousands of bytes composing digital images, with similar resource constraints.

From the many issues concerning WSNs, the problem of robustness against packet loss requires special attention. Indeed, in [START_REF] Yeo | Robust distributed multi-view video compression for wireless camera networks[END_REF], authors argue that packet losses are inherent in WSNs. Moreover, some real-world environments report loss rates surpassing the 70% of the transmitted packets [START_REF] Duran-Faundez | Error resilient image communication with chaotic pixel interleaving for wireless camera sensors[END_REF], thus error resilience is a critical issue.

Packet loss can occur due to many factors including channel errors, collisions and congestion. In WVSNs, these events, together with the vast amount of generated data (pixels), present problems to the communication process. Indeed, while to codify and communicate a single scalar measurement (such as temperature or presure) 2 or 3 bytes are largely enough, images are codified in thousand of bytes, grouped in hundreds of packets, which in case of getting lost adversely affect the quality of received images at the decoder side (more important, when image compression is applied), and their retransmission is power expensive.

Due to these statements, image transfer on WSNs requires communication techniques that maintain the quality of the obtained images, taking into account the problem of packet loss and traditional constraints of this kind of network. In some very constrained systems, traditional implementation of ACK or FEC-based methods is unsuitable because they require additional energy and time consumptions. Thus, it is important to evaluate alternative methods and adaptations of traditional techniques allowing more efficiently to deal with this issue.

A low-cost and effective solution for enhancing robustness of an image communication system is block interleaving. A block interleaving technique disrupts the traditional sequence of packetization-communication of an image blocks. It is based on the inherent spatial information redundancy of digital images. As neighboring blocks of pixels are sent in different data packets, there is a better probability of receiving neighboring blocks of a lost one, thus to count with enough information to estimate its value. In the literature, many works have adopted interleaving methods to increase robustness of a communication system but non of them compares these methods in a WVSNs context. Also, interleaving effectiveness seems to be assumed, and the evaluation of the robustness provided by the interleaving method is usually not rigorous. This paper addresses the problem of robustness in WVSNs evaluating various image interleaving methods found at the literature, addressing theoretical issues, and comparing them through extensive simulation, in terms of quality of the received images, with different testing parameters. The rest of this paper is organized as follows: Section II presents related works to the adoption and evaluation of image interleaving for robust image communication on WVSNs. Section III gives theoretical background on block interleaving, describing the selected block interleaving techniques studied in this work. Further analysis and parameters of evaluation for each technique are presented in Section IV. Simulation results are shown in Section V, presenting results for each method in terms of image quality for best parameter combinations found and results distributions associated to different result ranges. Finally, Section VI concludes this work.

II. RELATED WORKS

Various techniques can be adopted to allow reliability in a communication system. ARQ-based schemes can be adopted in order to ensure the reception of the packets sent by a camera source to the sink, at each hop in a path but this could be excessively expensive in terms of energy consumptions and delay. Energy savings can be achieved if a semi-reliable strategy is adopted [START_REF] Lecuire | Energy-efficient transmission of wavelet-based images in wireless sensor networks[END_REF], but the number of retransmissions needed to ensure the arrival of high priority packets can lead to similar consequences of a full-reliable communication, if no other error resilient techniques is applied. Similar consequences may occur by the transmission of redundant packets. For example, in [START_REF] Wu | Error resilient image transport in wireless sensor networks[END_REF], the problem of error probability in image transport is addressed by generating diverse copies of the same packet and transmitting them through diverse nodes in a cluster-based network. Copies are sent by a cluster head to various nodes, then, the copies are transmitted to the next cluster head towards the sink. Of course, if packets are lost, other copies could be received by the next cluster head. Such a technique ensures the reception of a good quantity of information packets at the decoder side, but the transmission of redundant packets and the required protocol complexity affects resource consumptions.

Multiple properties of digital images could be taken in advantage in order to provide resistant transport to communication errors in accordance with WVSNs specifications. In general, robust mechanisms proposed in the literature are based in the fact that, to many WVSNs applications, what is relevant is to obtain useful versions of the original images instead of lossless information. For example, the image communication scheme proposed in [START_REF] Wu | Multiple bitstream image transmission over wireless sensor networks[END_REF] takes advantage of this correlation. It is based on SPIHT coding of independent data blocks generated from parent-child relationships of wavelet coefficients. Each codified block generates an independent data stream, thus channel errors affect only the involved block of coefficients, reducing the impact of data loss.

In various works, inherent spatial correlation between blocks of pixels and their neighboring ones is exploited. In this paper, we address block interleaving as a technique for enhancing robustness of an image communication system. In the WVSNs context, work in [START_REF] El-Bendary | An efficient chaotic interleaver for image transmission over ieee 802.15.4 zigbee network[END_REF] shows a chaotic interleaver based on a discretized Baker map. New coordinates for a block are calculated by an equation taking as parameters original indices and a secret key vector [n 1 , . . . , n k ] with n 1 + • • • + n k = N , N being the number of blocks in an image row). This related work is interesting because it is specially conceived for Zigbee-based WVSNs, however, the enormous quantity of possible parameters combinations (secret key) makes difficult its exhaustive evaluation, thus this technique was not considered in the scope of this paper. Future studies should consider this technique.

In [START_REF] Duran-Faundez | Error resilient image communication with chaotic pixel interleaving for wireless camera sensors[END_REF], a chaotic Torus Automorphisms is adopted to improve robustness against packet loss in a Cyclops/Mica2 network. Results show that improvements in the quality of the reconstructed images can be achieved in presence of high packet loss rates, without excessively increase resource consumptions. Same technique is applied in [START_REF] Duran-Faundez | Tiny block-size coding for energy-efficient image compression and communication in wireless camera sensor networks[END_REF] for the case of images compressed in independent blocks of pixels. In this work we will evaluate three representative interleaving technique found in the literature.

III. BLOCK INTERLEAVING TECHNIQUES

Let us illustrate the effect of packet loss in the transmission of a digital image with the example in Figure 1.

The captured image of Figure 1a is received with losses, obtaining the bitmap of Figure 1b. After applying an error concealment method (as those mentioned by [START_REF] Turner | Image transfer: an end-to-end design[END_REF]) it is possible to obtain the reconstructed version of Figure 1c. Notably, the quality of the received image is greatly degraded. In cases where lost data packets are isolated a good reconstruction can be achieved (see, e.g., the half-bottom of the reconstructed image), but in places where big areas are lost (in some cases due to burst losses), reconstruction is impossible without appropriate codification at the source side. Block interleaving can be applied to mitigate packet loss effects. The principle of block interleaving is to transfer adjacent blocks of pixels in the original image (codified or not) into data packets as far as possible. The ideal situation is not finding two neighboring blocks from the original image in the same packet, so there is better probability of receiving neighboring blocks to a lost one (and thus, of better estimate the missing values). Virtually, an interleaving algorithm performs this task disrupting the bitmap of the original image, assigning different positions to adjacent pixels, and thus increasing the probability of transmitting in different packages, as illustrated in Figure 2.

For a similar scenario of the illustrated in Figure 1, the image in Figure 2a is sent through a non-reliable channel, but this time being interleaved first. The interleaved version of the original image is shown in Figure 2b. At the decoder side, interleaved bitmap suffers same packet loss of previous example, as shown in Figure 2c. After performing inverse interleaving (reallocation of blocks of pixels to the original position), we obtain the bitmap of Figure 2d. In this new map, missing information is noticeably better isolated, giving as a result (after the error concealment application), the approximate version of Figure 2e. This example illustrates the improvement allowed by a block interleaving in comparison with the sequential transmission case.

Of course, the obtained results will depend on the adopted interleaving technique. In the following, we explain and evaluate three interleaving methods applied in the literature for robust image communication.

A. Turner and Peterson (TP) interleaving

The encoding algorithm introduced by Turner and Peterson (TP) [START_REF] Turner | Image transfer: an end-to-end design[END_REF] consists into modify the original spatial positions of each pixel by two parameters. The first, byteOffset, specifies the distance between adjacent pixels (in the original image) into a single packet, i.e., it specifies the number of pixels not included between two pixels added to a packet. The second parameter, packetOffset, indicates the distance (in pixels) between the first elements of two adjacent packets, starting the generation of a new one packetOffset pixels after the first element of the previous packet, diminishing the effects of burst packet loss. In the paper, this method is not applied to WSNs, but given its simplicity, its implementation is perfectly possible to this kind of network. Then, for I(H × W ) images, I = {p 0 , p 1 , . . . , p H * W -1 }, where p corresponds to the pixels of I; the i-th packet P i , such that i ≥ 0, contains m elements, being Pi = {pj , . . . , pj+2 * ByteOffset, . . . , p j+(m-1) * ByteOffset }, with j = i * PacketOffset.

B. DSJ-AL interleaving method

As a part of the Dual Stream JPEG Adaptation Layer, DSJ-AL, method, Posnak et al. [START_REF] Posnak | Techniques for resilient transmission of jpeg video streams[END_REF] incorporate an interleaving method in which, for an image of H × W blocks, two originally adjacent blocks are separated by a parameter step. This is performed as follows: Given a packetized block in coordinate (x i , y i ) (it is assumed that the first position is (0, 0)), the coordinate of the following block to be packetized, (x i+1 , y i+1 ), is calculated by using Equations 1 and 2.

y i+1 = (y i + step) mod W (1)
x i+1 = x i + y i + step W (2)
Originally, this technique is applied to interleave codified blocks of pixels, but it can also be used to interleave individual non-codified blocks. For proper generation of new positions, it is indicated that, once the maximum number of rows (H) for x i+1 the next pixel to be assigned will be the following to the first on the current cycle. For example, assuming that the first coordinate in the current cycle is (0, 0), once x i+1 > H exceeds this amount, the coordinate generated for the next pixel should be (0, 1).

C. Torus Automorphims (TA) interleaving

Torus Automorphims (TA) are highly chaotic systems, that can be used as permutation functions in order to implement block interleaving. In this work, we evaluate Arnold's Cat map-based automorphisms applied in [START_REF] Duran-Faundez | Error resilient image communication with chaotic pixel interleaving for wireless camera sensors[END_REF]. Given an image I of N × N blocks, and {x, y ∈ N / 0 ≤ x < N ∧ 0 ≤ y < N }, each pixel with coordinates (x, y) is stored in a new position (x ′ , y ′ ) calculated by Equation 3.

x ′ y ′ = 1 1 k k + 1 n x y mod N (3)

IV. FINDING EVALUATION PARAMETERS

An important topic in this work is to evaluate the performance of the selected interleaving methods. In image communication through WVSNs these techniques are used as a tactic to obtain robustness against losses, thus their performance is directly related to the quality of the resulting images at the receiver side. Specifically, the stated problem consists on the choice of one (or various) interleaving techniques and its parameters, so that the quality of the obtained images is within acceptable ranges. Due to the above, we study the behavior of the described techniques with multiple images and loss patterns, using the simulation framework Sim-LIT [START_REF] Orellana-Romero | Sim-lit: A simulation framework for image quality assessment in wireless visual sensor networks under packet loss conditions[END_REF]. To assess conditions close to reality, 122 loss patterns obtained from real-world experiments (see full description in [START_REF] Duran-Faundez | Error resilient image communication with chaotic pixel interleaving for wireless camera sensors[END_REF]) and monochromatic (8bpp) 128×128 versions of the 44 test images found at [START_REF] Sipi | Sipi image database -volume 3: Miscellaneous[END_REF] are used. The following subsections describe the adopted range of parameters used in our simulations.

A. TP method parameters

In TP, there is a large number of parameter combinations available for selection, making difficult the choice of good combinations. Moreover, when testing TP with arbitrary parameters, for some particular parameter combinations, coordinates generated by the algorithm are equal for various different pixels, overwriting the first ones.

In order to obtain a list of valid parameter combinations, we observe the reallocation of pixels in a 128 × 128 boolean matrix, initiated to False. Then, generating coordinates with TP technique, True values are assigned to each new generated position, labeling the combination of parameters as invalid if any overwriting occurs. This study is performed for combinations of parameters included in the ranges 0 ≤ ByteOffset < 16, 500 and 0 ≤ PacketOffset < 16500, using packet size of 27 bytes. Finished this evaluation, the amount of correct pairs is greater than 8300, being necessary to check the existence of identical interleaving patterns generated using different combinations of parameters. After eliminating equivalent combinations (parameters who through same interleaving patterns that other), valid pairs are decreased to 8192. We must note that pairs qualified as valid parameters are bounded by a value of 128 × 128, equivalent to the length of adopted images. From the previous study, we obtained a table with valid pairs within the specified range.

B. DSJ-AL interleaving parameters

In order to evaluate the behavior of DSJ-AL interleaving, we use a 128 × 128 image. Simulations varying its unique parameter step in the range [START_REF] Akyildiz | Wireless sensor networks: A survey[END_REF]17000] were performed, calculating the PSNR between the mixed and the original images. Infinity results (cases where resulting and original images are equal) are obtained for step = 1 and step ≥ 128 × 128.

Comparing the results with DSJ-AL behavior, we note coherence between them, because, for step = 1 there is no offset in the pixels, and for step ≥ 128 × 128 a greater distance between pixels to the total length of the image, is generated. This behavior generates unmixed images for step = 1 and for step ≥ (N × N ) images. In this way, the evaluation range is defined to 1 < step < (N × N ), N ∈ N, for a N × N image.

C. TA interleaving parameters

Because of the presence of two parameters (k and n), there is a great amount of combinations, making difficult to determine a parameter pair which generates good results. One aspect to consider is that the technique is based on toroids. This suggests that the results are generated in cycles, allowing narrow each parameter to a range. Generally speaking, it is expected that there is a maximum value k with its corresponding n value, that limits the results, reducing the amount of possible choices. Thus, the method used to define limits is based on Equation 3, where it is possible to note that the interleaved positions (x ′ , y ′ ) are determined by the matrix 1 1 k k + 1 n mod N . Indeed, when this matrix corresponds to the identity, it causes x ′ = x and y ′ = y, generating the same original image after interleaving. To determine the parameter values that enforce this condition, Equation 3 is evaluated for an arbitrary range of values (not considering n = 0), storing the first parameter n = T that generates the identity matrix, for each value of k in a bounded parameter list. This approach is carried out for images of size 128 × 128 pixels, evaluating the parameters for the range 0 ≤ k ≤ 300 y 1 ≤ n ≤ 300.

When analyzing the generated data, a pattern for T values when varying k is determined defined by T (k) = T (k + R × 128), with R ∈ N y 0 ≤ k ≤ 128. That is, the sequence of values of n that generate an original image for different values of k is repeated each 128 increments from this last one. Given this, it is possible to establish as limit for k the value 128 in images of 128 × 128 pixels. For the parameter n, it is detected that the values for which we obtain the same original images n = T are discrete [START_REF] Soro | A survey of visual sensor networks[END_REF][START_REF] Yeo | Robust distributed multi-view video compression for wireless camera networks[END_REF][START_REF] Lecuire | Energy-efficient transmission of wavelet-based images in wireless sensor networks[END_REF][START_REF] Posnak | Techniques for resilient transmission of jpeg video streams[END_REF][START_REF] Wu | Multiple bitstream image transmission over wireless sensor networks[END_REF]16,24,32,48,64,96,128) and they have as superior limit 128, so that, for each value of k, in images of 128 × 128 pixels, an original bitmap is obtained, at least once, when n = 128.

After performing similar experiments with N = 32, N = 64 and N = 256, it is possible to affirm that for square images of N × N blocks, with N = 2 x / x ∈ N, the range of k values producing different interleaved patterns can be limited to 0 ≤ k ≤ N . Moreover, it is observed that the discretized nature of T is maintained, adding, or eliminating, values to the set of frequencies, for a given N . In this way, for N = 32, resulting T values are in the set F 32 = {3, 4, 6, 8, 12, 16, 24, 32}, for N = 64, F 64 = F 32 {48, 64}, for N = 128, F 128 = F 64 {96, 128} and for N = 256 they are given by the set F 256 = F 128 {192, 256}. Generalizing, to image interleavings with N = 2 x / x ∈ N, the set of frequencies for value T are limited by T = N and corresponds to the set

F N = F N =2 x-1 {2 x-1 + 2 x-2 , N } / x > 2.

V. SIMULATION RESULTS

Once the range of different parameter combinations for evaluation is defined, we proceed to assess which configuration(s) allows generate better results when recovering a transmitted image, so we can count on the needed information to obtain the required results. For this, the following procedure is performed for each method, evaluating obtained results looking for the best parameters in each case.

A. Generic evaluation process

Using the described set of 44 images and 122 loss patterns, we simulate with Sim-LIT the transfer of each of the images for each of the available loss patterns, using for that interleaving methods with parameters combinations as described in Section IV (for H = W = N = 128), depending on the employed method. Average value of the correctly received neighboring pixels to a lost one was adopted as error concealment technique. between the resulting and the original images are obtained, for each simulation, assessing image quality for the different packet loss patterns and for the different configurations of the interleaving methods. We define a packet size of 27 bytes which corresponds to the default payload value of the Cyclops camera firmware. Once obtained these data, parameters providing the best results are evaluated. To this, the mean PSNR of the 44 reconstructed images is calculated for each simulation. In this way, a set of PSNR measures is obtained to measure the performance of each parameter combination under different channel conditions (loss patterns). Once the set of values that represents the different valid combinations of parameters is generated, it is necessary to compare their performance to determine a combination who gives the best results. To this, we adopted GNU Octave software, generating a curve representing the obtained points after simulation, adopting to that polynomial functions of degree 8, with function polyfit(). Thus, to determine the best pair of parameters the area under each curve is calculated, using functions polyout() to construct the polynomial function (using calculated coefficients with polyfit()), and quad() to calculate the defined integral in the range 1 to 90 (loss rates), determining the biggest area under each curve, and thus, the parameters generating the best results. We must note that this evaluation technique is not valid in every case, because it adopts an approximate representation method. However, the obtained results from this analysis are valid for most of the studied cases. To analyze the data produced after simulation in a generic manner, interpolated curves for each valid parameter are generated and the area under each of these is calculated. After that, the range of obtained areas is divided in 30 intervals and histograms are constructed with these data, showing the number of parameters for which areas are obtained within each group.

B. Parameter evaluation of Torus Automorphims

After performing the described procedure, the parameters providing a better value when calculating the integral of the interpolated curves are k = 28 and n = 114, obtaining an area under the curve of 2681. In the same way, parameters k = 126 and n = 2 generate the minimum area value being of 2194. In the graphics of Figure 3, simulation data, along with interleaved curves for the best and worst parameters combinations is shown.

C. Parameter evaluation for DSJ-AL interleaving

DSJ-AL interleaving has a unique parameter, s. Following the procedure of Section V-A, it is obtained s = 92 as the best parameter, reaching a maximal area under the interpolated curve equal to 2657. Moreover, the minimal area is determined, reached by the interpolated curve calculated with parameter s = 15082, with a value of 2116. Graphics of Figure 4 shows simulation data, along with interleaved curves for the best and worst parameters combinations.

D. Parameter evaluation for TP method

TP method counts on two parameters, ByteOf f set y P acketOf f set. After applying the evaluation procedure, the pair giving the best results, measured through the calculation Figure 5 shows the respective results.

E. Results comparison

To compare the evaluated interleaving methods and a sequential image transmission (without interleaving), Figure 6 shows all obtained data using the best parameters of each method and the data generated when a sequential (noninterleaved) communication is applied. From this, it is possible to say that the all evaluated interleaving techniques generate better average results than the sequential transmission, when packet loss rates surpass 15%. Below this value, all transmission strategies are quite similar. Regarding interleaving techniques, it is clear that generate similar results, however, there are differences over 40% of packet loss. Thus, for loss rates in the range between 40% and 80% TA generates slightly better results than their peers, however, TP method throw quite near results. Meanwhile, DSJ-AL is the one with a lower performance compared to other interleaving methods. Table I summarizes the maximum, minimum, and mean of the areas under the interpolated curves.

Notably, these are the maximum for each of the technical values, so to generate a conclusion on the best performance is necessary to evaluate the above results as a whole. Histograms of Figure 7 show the different distributions of the quantities of parameters generating areas under the interpolated curves for the different defined intervals. From the results, we observe that the methods having a bigger quantity of parameters giving results near the maximum are TA and TP. In the case of DSJ-AL, the major quantity of results are near the minimum.

Comparing the best evaluated methods, we can say that TA stands out above TP, because its distribution throws a larger quantity of parameters giving areas near the maximum.

VI. CONCLUSION

The conclusions generated are defined based on image transfer simulations for different cases, considering the quality of the final image as evaluation metric, measured by PSNR. Analyzing the simulated results for each technique, their methods of operation, and better operating range parameters are determined in terms of the quality of the final image in the receiving node. As a result, we observe that all these methods give better results than the sequential transmission for a loss rate greater than 15%. Below this value the results of the final image qualities are very similar for all the evaluated cases. Regarding the difference between the interleaving algorithms, it is determined that gaps are generated in the final quality of the images transmitted with a packet loss superior to 40%. For the evaluation of TA is determined that the parameters that generate better quality in the final image correspond to k = 28 and n = 114, while for evaluation Turner and Peterson interleaving is determined that the parameters that generate better quality in the final image are byteOf f set = 11727 and packetOf f eset = 5333, considering as indicator the integral of the curve interpolated by a polynomial of degree 8 for each of the simulated data. For the evaluation of DSJ-AL interleaving results show that the best performing parameter in terms of the quality of the final image is s = 92, based on the same study generated for the above methods. Generally, after comparing the number of parameter combinations that generate near maximum delivered by each best result for the methods studied and compare the interpolated curves generated for each combination of them better results, it is determined that TA interleaving maintains the best performance images for different channel states, being nearly followed by TP interleaving.

Fig. 1 .

 1 Fig. 1. Example of the effect of packet loss on sequential image transmission.

Fig. 2 .

 2 Fig. 2. Example of the application of image interleaving in an image communication scenario in presence of high packet loss.

Fig. 3 .

 3 Fig.3. Results of the best and worst pair of parameters, using TA interleaving.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Results of the best and worst pair of parameters, using DSJ-AL interleaving.

Fig. 6 .

 6 Fig.6. Results for sequential and interleaved image transmission using best parameters.

Fig. 7 .

 7 Fig. 7. Histogram of areas under the curve for each method.

ACKNOWLEDGMENT

This work is supported by Conicyt through Fondecyt Project No. 11121657.


J Lambert

H Chouh

G Rougeron

V Bergeaud

S Chatillon

L Lacassagne

J C Iehl

J P Farrugia

V Ostromoukhov

Interactive Ultrasonic Field Simulation For Non-Destructive Testing

Keywords: Descriptors (according to ACM CCS): F, 1, 2 [Computation by Abstract Devices]: Modes of Computation-Interactive and reactive computation

In Non Destructive Testing domain there a is a need for interactive simulation tools. This poster presents such a tool which computes up to 20 ultrasonic field images per seconds for simple configurations on a high-end CPU. Extensions to handle more complex configurations are also presented.

Introduction

Inspection simulation is used in a lot of non-destructive testing (NDT) applications in industry. The aim of Civa Software is to provide fast and easy to interpret NDT simulations and reconstructions. In this context, this poster presents a new tool for the interactive simulation of ultrasonic fields.

Ultrasonic Field Simulation Model

An ultrasonic field simulation computes the field radiated by a probe (immersed or in contact, mono-element or phased array) in an inspected specimen. It is generally used to choose the best suited characteristics of a probe for a given control. Ultrasonic waves propagation into solids involves complex physical phenomena. Among them, the following ones can be highlighted. Different types of volumic waves can propagate at different celerities. At each surfacic interface they can be specularly reflected or refracted, with possibly wave type conversion. Specimen can be homegeneous or inhomogeneous. Materials can be isotropic or anisotropic. † In Civa Software, an ultrasonic field simulation consists in three steps [Gen03]:

1. performing pencil tracing by computing for each pencil on one hand an axial ray that starts from a source point † Anisotropy implies the waves celerity depends on the propagation direction.

on the probe and reaches the discretized field area, and on the other hand the following parameters : amplitude (via a divergence factor and the multiplication of Fresnel coefficients), time of flight and duration. 2. delaying and accumulating pencil contributions on a vectorial displacement impulse response :

---→ IR D (t) ∈ C 3 according to the delay laws used for the probe, 3. convoluting this response with the probe signal and extracting a max amplitude for each field point which leads to a colored image.

Field simulations in current Civa software can handle very complex configurations at the price of often important computation time (typically from dozens of seconds to hours). A new tool adressing simple configurations has therefore been developped to reach our interactivity goal.

A fast implementation for simple configurations

A simplified model, restricted to homogeneous isotropic specimens with planar boundaries has been developped. Single refracted ray paths are computed through a 1D-Newton equation resolution ( [WZHB09]). A supplementary reflection without conversion can easily be added. Pencil parameters are computed directly through analytic formulae.

A first reference, OpenMP parallelized, scalar, single precision implementation in three steps, using Intel MKL FFT, has been coded on CPU. It has then been further accelerated. Using Boost.SIMD library, steps 1. and 3. have been submitted to EUROGRAPHICS 2014. vectorized. Morevover, a single step algorithm avoiding the unnecessary storage of temporary data and providing more work for each thread has been settled. The resulting optimized code scales well on a 2x12 cores CPU. Its overall performance is around 3.5x faster than previous reference implementation on a set of representative configurations. For the simplest ones, it reaches 20 fps. This implementation shall be tested on other harwdares (CPUs with AVX/AVX2 registers and Intel MIC). An optimized GPU Cuda based code shall be implemented.

Towards more complexity

Previous implementation can be functionnaly extended in order to take into account :

• non planar surfaces, which raises the complexity of nonlinear equations to solve, • the length of ray paths with a second reflection or refraction on the specimen boundaries with possibly wave conversion, which leads to solve a system of equations • specimen with many boundaries, which means some heuristic must be used in order to avoid testing all possible combinations of surfaces for a given path type.

[JM12] might be used to find neighbor paths starting from a known specular path.

But as these extensions remains limited, for even more complex configurations, i.e. :

• specimen with any geometry (non-planar surfaces or 3D triangle mesh) • inhomogeneous specimen with many different materials • anisotropic materials • paths of any length another method based on fast ray tracing has to be developped. Early experimentations showed that an ultrasonic raytracer based on Intel Embree on CPU or Nvida OptiX on GPU should allow us to reach the necessary level of performance required by interactivity (at least via progressive computation). Furthermore a geometrical method has been settled to trace rays in anistropic materials much faster than Civa analytic method. Fast methods for Fresnel coefficient and divergence factor computation might be worked out in order to compute full pencil parameters.

Relying on these fast tools, the main algorithmic effort shall concern the pencil computation step as for complex configurations it represents by far the major part of computation time. Two different ideas shall be explored :

1. Pencil computation in two sub-steps :

• A geometrical sub-step = For each field point and each path type, find solid angles through which the probe surface is seen via a reverse ray tracer guided by an adaptive algorithm (an optimization algorithm aiming at minimizing the distance of rays to the probe surface, or AD-Frustum [CLT * 08]) • A physical sub-step = For each field point, each path type and each solid angle, discretize the surface probe into a set of source points and perform pencil tracing starting from the probe.

2. Pencil computation in one step : Perform a MonteCarlo light tracing (or Markov Chain Monte Carlo) with importance sampling, and possibly regularization (like in [KD13]), following ray paths starting from the probe until they reach the field area and contribute to the field points impulse responses.

Perspective and conclusion

A fast US field simulation tool allowing interactive change of any parameter for simple configurations has been developped on CPU. New hardware shall be tested (CPUs with larger SIMD registers, Intel MIC, and GPUs). The method shall be extended to take into account more complex geometries or longer ray paths. For very complex configurations, in particular those with anisotropic materials, an alternative method using ray tracing is currently being developped. The resulting interactive tools will be integrated in next Civa software commercial release. As different codes developed for different hardware will finally be available, an auto-tuning mechanism will be settled in order to choose automatically the best one for a given ultrasonic field configuration. Another mechanism might be developed to automatically tune computation options in order to keep a satisfactory level of interactivity [GD12]


Taro Kimura
email: taro.kimura@cea.fr

Duality and integrability of a supermatrix model with an external source

We study the Hermitian supermatrix model involving an external source. We derive the determinantal formula for the supermatrix partition function, and also for the expectation value of the characteristic polynomial ratio, which yields the duality between the characteristic polynomial and the external source with an arbitrary matrix potential function. We also show that the supermatrix integral satisfies the 1-and 2D Toda lattice equations as well as the ordinary matrix model.

. . .

Introduction

The matrix model has played a crucial role in quantum field theory (QFT) for the last few decades, as a 0D QFT model, or rather a toy model for the infinite-dimensional path integral. Since it is just defined by the finite-dimensional integral of a matrix itself, one can obtain an exact solution in various situations, which provides a significant insight for the understanding of QFT. This methodology is also referred to as random matrix theory (RMT), and is now applied to an extremely wide range of research fields [START_REF]The Oxford Handbook of Random Matrix Theory[END_REF].

In QFT, in order to compute correlation functions, it is convenient to introduce the generating function by adding the extra source term Z[J ] = Dφ e -1 S[φ]+ d D x J(x)φ (x) .

(

) 1
The correlation function can be obtained by taking the functional derivative with the source field

φ(x 1 ) • • • φ(x k ) = δ δ J (x 1 ) • • • δ δ J (x k ) log Z[J ] J =0 . ( 2
)
This generating function is defined in the sense of a path integral, and thus it is quite difficult to compute in general. On the other hand, the matrix model with an external source, which is just given by a finite-dimensional integral, plays a similar role to the generating function in QFT,

Z N (A) = d X e -1 TrW (X )+TrX A , ( 3
)
where the integral is taken over a size-N matrix. We can compute the correlation function by taking the derivative with respect to the source matrix A, as well as the QFT generating function.

PTEP 2014, 123A01

T. Kimura A remarkable feature of the matrix model with an external source is the duality with a characteristic polynomial, which was found in Ref. [2] especially for the Gaussian matrix model. The claim of the duality is as follows: The M-point correlation function of characteristic polynomials in the size-N matrix model with an external source is equivalent to the N -point function with a size-M matrix, by exchanging the arguments of the characteristic polynomial and the external source. This duality is now extended to the generic β ensemble [3], the chiral ensemble [4], and beyond the Gaussian model [5], and has various interesting interpretations in terms of conformal field theory, string/ M-theory, knot theory, algebraic geometry, and so on.

Another interesting property of the external source matrix model is the connection with the integrable system. The matrix integral with an external source and also the characteristic polynomial average satisfy the integrable equation, e.g., the Toda lattice equation, and they can be interpreted as the corresponding τ -function; see, e.g., Ref. [6]. Such integrability of the matrix model is quite useful to determine the spectral density and also the correlation function.

In this paper we study a supersymmetric version of the matrix model involving an external source based on a Hermitian supermatrix. The supermatrix method is now well known both in high-energy and condensed-matter physics. For example, one can avoid technical difficulty in dealing with the normalization factor of the correlation function by applying the supermatrix method, in a similar way to the replica method. This property helps us to take the disorder average in a random system [START_REF] Efetov | Supersymmetry in Disorder and Chaos[END_REF]. In Sect. 2 we will compute the U(N |M) symmetric supermatrix partition function involving the external source. As in the case of the ordinary Hermitian matrix model, we can utilize the determinantal structure in order to perform the integral. We will obtain a determinantal formula for the partition function, which can be expressed as a size-N determinant of the U(1|1) partition function as a determinantal kernel especially for N = M.

As in the case of the ordinary matrix model, there have been some proposals for the duality of the supermatrix model. For example, in particular for the Gaussian supermatrix model, one can show an explicit duality between the external source and the characteristic polynomial [START_REF] Desrosiers | [END_REF]. In Sect. 3 we will exhibit that this duality generally holds for the supermatrix with an arbitrary potential function, and is just interpreted as a Fourier transformation. In addition, it will be pointed out that the relation between the characteristic polynomial ratio with the ordinary matrix model [9] and supermatrix models [10][11][12] is naturally explained as some specialization of the characteristic polynomial expectation value in the supermatrix.

In Sect. 4 we will also show the integrable equations for the supermatrix integral with the external source and the characteristic polynomial, based on the determinantal formula derived in this article. Since we have more parameters for the partition function than the ordinary matrix model, we can obtain several integrable equations corresponding to variables parametrizing the external source and the characteristic polynomial.

Supermatrix model with external source

Let us introduce a Hermitian supermatrix model with a potential function W (x), involving the external source matrix C, as a natural generalization of the Hermitian matrix model (3),

Z N ,M {a i } N i=1 , {b j } M j=1 = d Z e -1 StrW (Z )+StrZC , (4)
PTEP 2014, 123A01

T. Kimura where is the coupling constant. The Hermitian supermatrices Z and C are given by

Z = X ξ ξ † Y , C = A η η † B , ( 5
)
where X , A are N × N , and Y , B are M × M (bosonic) Hermitian matrices. ξ and η are N × M fermionic matrices, whose elements are given by Graßmannian variables. We can assume that C is a diagonal matrix without loss of generality, and in this paper we also assume N ≥ M. The matrix measure of ( 4) is invariant under the supergroup transformation, d Z = d(U ZU -1 ) with U ∈ U(N |M), and is normalized by the corresponding volume factor Vol. U(N |M). As well as the ordinary Hermitian matrix integral, one can decompose this measure into diagonal and angular parts,

d Z = N ,M (x; y) 2 d N x d N y dU, ( 6
)
and the corresponding Jacobian N ,M (x; y) is now given by the generalized Cauchy determinant [13],

N ,M (x; y) ≡ N (x) M (y) N i=1 M j=1 (x i -y j ) -1 = det x k-1 i (x i -y j ) -1 with ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ i = 1, . . . , N j = 1, . . . , M k = 1, . . . , N -M , ( 7
)
where N (x) is the Vandermonde determinant,

N (x) = det 1≤i, j≤N x j-1 i = N i< j (x i -x j ). ( 8
)
The fact that the Jacobian factor can be written as a determinant will play a crucial role in the following discussion.

To perform the angular part integral of (4), we now apply the supergroup version of the Harish-Chandra-Itzykson-Zuber (HCIZ) formula, which is written in terms of the Cauchy determinant [14][15][16],

I N ,M (Z , C) ≡ U(N |M) dU e StrZUCU -1 = 1 N ,M (x; y) N ,M (a; b) det 1≤i, j≤N e x i a j det 1≤i, j≤M e -y i b j . ( 9
)
This is a quite natural extension of the original formula for the U(N ) group, which is given by replacing the Cauchy determinant with the Vandermonde determinant [17,18]:

I N (X, A) ≡ U(N ) dU e TrXU AU -1 = 1 N (x) N (a) det 1≤i, j≤N e x i a j . ( 10
)
Applying the formula (9), we obtain the expression only in terms of eigenvalues:

Z N ,M = 1 N ,M (a; b) N i=1 dx i e -1 W (x i )+x i a i M j=1 dy j e 1 W (y j )-y j b j N ,M (x; y). ( 11
)
PTEP 2014, 123A01

T. Kimura

Let us comment on the convergence of this integral. Naively speaking, it is not a convergent integral in general. We can regard this as a formal integral, but one possible way to avoid this divergence is to modify the matrix potential by inserting a constant supermatrix

I N ,M = diag(1 N , i1 M ) ∈ U(N |M), i.e., StrZ 2 = TrX 2 -TrY 2 → StrI N ,M Z 2 = TrX 2 + TrY 2 [3]
. However, this treatment does not work for the generic potential function W (x). Another possibility is to consider the coupling constant as an imaginary number, ∈ iR. In this case, although the integrand becomes an oscillating function, one can provide an interpretation of it as the Fresnel integral. For example, this interpretation works well for the Chern-Simons matrix model [19] and its supermatrix generalization, which is called the ABJM matrix model [20][21][22], although the corresponding coupling constant is pure imaginary, e.g.,

= 2πi/(k + N ), or = 2πi/k.
Then, due to the Cauchy determinantal formula [START_REF] Efetov | Supersymmetry in Disorder and Chaos[END_REF], it can be written as a size-N determinant, consisting of N × M and N × (N -M) blocks,

Z N ,M {a i } N i=1 , {b j } M j=1 = 1 N ,M (a; b) det Q k-1 (a i ) R(a i ; b j ) , ( 12
)
where the indices run as i = 1, . . . , N , j = 1, . . . , M, k = 1, . . . , N -M, and we have introduced the following functions:

P i-1 (x) = x i-1 , R(x; y) = 1 x -y , ( 13
)
Q i-1 (a) = dx P i-1 (x) e -1 W (x)+xa , (14)
R(a; b) = dxdy R(x; y) e -1 W (x)+ 1 W (y)+xa-yb . ( 15
)
Note that the function (14) obeys

Q i (a) = d da Q i-1 (a), ( 16
)
and Q i=0 (a) is seen as a generalized Airy function:

Q i=0 (a) = dx e -1 W (x)+xa . ( 17
)
As we will see later, these functions are just seen as a (double) Fourier or rather Laplace transform of P i-1 (x)e -1 W (x) and R(x; y)e -1 W (x)+ 1 W (y) , respectively. 1 Let us comment on some specialization of the formula ( 12). If we take the limit M = 0, the determinantal expression ( 12) is reduced to the the well known formula for the ordinary Hermitian matrix 1 If we modify the source term in the supermatrix integral

Z N ,M {a i } N i=1 , {b j } M j=1 = d Z e -1 StrW (Z )+i StrZC , ( 18
)
they just become Fourier transforms:

Q i-1 (a) = dx P i-1 (x) e -1 W (x)+i xa , R(a; b) = dxdy R(x; y) e -1 W (x)+ 1 W (y)+i xa-iyb . ( 19
) 4/12
PTEP 2014, 123A01

T. Kimura model with the external source

Z N {a i } N i=1 = d X e -1 TrW (X )+TrX A = 1 N (a) det 1≤i, j≤N Q j-1 (a i ); (20)
see, e.g., Ref. [6] for details. Then, in the case with M = N , the partition function becomes

Z N ,N {a i } N i=1 , {b j } N j=1 = 1 N ,N (a; b) det 1≤i, j≤N R(a i ; b j ). ( 21
)
This shows that the formula for U(N |N ) theory is reduced to a size-N determinant of the "dual" kernel function R(a; b), which is given by the U(1|1) partition function

R(a; b) = Z 1,1 (a, b) a -b . ( 22
)
This factorization property is called Giambelli compatibility [23], and follows the Fay identity [START_REF] Fay | Theta Functions on Riemann Surfaces: Lecture Notes in Mathematics[END_REF].

Characteristic polynomial average and duality

In this section we consider expectation values of the characteristic polynomial ratio with the supermatrix model in the presence of an external source:

N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 = d Z e -1 StrW (Z )+StrZC p α=1 Sdet (λ α -Z ) q β=1 Sdet μ β -Z -1 , ( 23
)
and then show the duality between the ( p|q)-point function with U(N |M) theory and (N |M)-point function with U( p|q) theory. The characteristic polynomial average (23) includes several useful situations:

• M = q = 0: Characteristic polynomial with the Hermitian matrix model

N ; p {a i } N i=1 ; {λ α } p α=1 = d Xe -1 TrW (X )+TrX A p α=1 det(λ α -X ). ( 24
)
• M = p = 0: Characteristic polynomial inverse with the Hermitian matrix model

N ; q {a i } N i=1 ; {μ β } q β=1 = d X e -1 TrW (X )+TrX A q β=1 det(μ β -X ) -1 . ( 25
)
• M = 0: Characteristic polynomial ratio with the Hermitian matrix model

N ; p,q {a i } N i=1 ; {λ α } p α=1 , {μ β } q β=1 = d X e -1 TrW (X )+TrX A p α=1 det(λ α -X ) q β=1 det(μ β -X ) -1 . ( 26
)
Therefore, expression (23) provides a master formula for the characteristic polynomial average in various matrix models. For example, the duality in the case with M = 0 (26) claims that the characteristic polynomial ratio with the ordinary Hermitian matrix model is dual to another supermatrix model, as shown in Ref. [11].

PTEP 2014, 123A01

T. Kimura

As discussed in Sect. 2, the angular part of the integral is performed using the HCIZ formula ( 9), and then we obtain the expression only in terms of eigenvalues:

N ,M; p,q = 1 N ,M (a; b) N i=1 dx i e -1 W (x i )+x i a i M j=1 dy j e 1 W (y j )-y j b j N ,M (x; y) × p α=1 N i=1 (λ α -x i ) M j=1 (λ α -y j ) q β=1 M j=1 (μ β -y j ) N i=1 (μ β -x i ) . ( 27
)
It is now convenient to apply the following identity to this expression:

N + p,M+q (x, λ; y, μ) = N ,M (x; y) p,q (λ; μ) p α=1 N i=1 (λ α -x i ) M j=1 (λ α -y j ) q β=1 M j=1 (μ β -y j ) N i=1 (μ β -x i ) . ( 28
)
If we assume N + p > M + q, the LHS of this equation is a large Cauchy determinant of the size N + p, with several block structures,

N +q,M+q (a, λ; b, μ) = det ⎛ ⎜ ⎝ x k-1 i λ k-1 α x i -y j -1 λ α -y j -1 x i -μ β -1 λ α -μ β -1 ⎞ ⎟ ⎠ (29) = (-1) N p+Mq det ⎛ ⎜ ⎝ λ k-1 α x k-1 i λ α -μ β -1 x i -μ β -1 λ α -y j -1 x i -y j -1 ⎞ ⎟ ⎠ , ( 30
)
where i = 1, . . . , N , j = 1, . . . , M, α = 1, . . . , p, β = 1, . . . , q, and k = 1, . . . , N + p -Mq.

From this expression, we obtain the determinantal formula for the characteristic polynomial ratio expectation value

N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 = 1 N ,M (a; b) p,q (λ; μ) N i=1 dx i e -1 W (x i )+x i a i M j=1 dy j e 1 W (y j )-y j b j N + p,M+q (a, λ; b, μ) = 1 N ,M (a; b) p,q (λ; μ) det ⎛ ⎜ ⎝ Q k-1 (a i ) P k-1 (λ α ) R(a i ; b j ) S R (λ α ; b j ) S L (a i ; μ β ) R(λ α ; μ β ) ⎞ ⎟ ⎠ , ( 31
)
where we have introduced the following auxiliary functions in addition to (13), (14), and (15):

S L (a; μ) = dx 1 x -μ e -1 W (x)+xa , S R (λ; b) = dy 1 λ -y e 1 W (y)-yb . ( 32
)
The formula (31) shows a duality between the external source and the characteristic polynomials in the supermatrix model with an arbitrary potential function W (x), since it can 6/12

PTEP 2014, 123A01

T. Kimura also be written as follows:

N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 = (-1) N p+Mq N ,M (a; b) p,q (λ; μ) det ⎛ ⎜ ⎝ P k-1 (λ α ) Q k-1 (a i ) R(λ α ; μ β ) S L (a i ; μ β ) S R (λ α ; b j ) R(a i ; b j ) ⎞ ⎟ ⎠ . ( 33
)
This duality was originally shown in the Gaussian model with the harmonic potential W (x) = 1 2 x 2 [START_REF] Desrosiers | [END_REF]. Under this duality, we see the correspondence between the auxiliary functions:

P k (λ) ←→ Q k (a), S L (a; μ) ←→ S R (λ; b), R(a; b) ←→ R(λ; μ), ( 34
)
and then the ( p|q)-point function for the U(N |M) supermatrix model with the potential W (x) is converted into the (N |M)-point function for the U( p|q) supermatrix model with another potential W (x), which is obtained through Fourier transformation:

N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 F.T. = p,q;N ,M {λ α } p α=1 , {μ β } q β=1 ; {a i } N i=1 , {b j } M j=1 . ( 35
)
In this sense, it is easier to see this duality, especially for the Gaussian matrix model, because the Gaussian function is self-dual with respect to Fourier transformation. We remark that this kind of duality can also be found in the Chern-Simons matrix model and the ABJM matrix model, which are obtained by replacing the Vandermonde determinant with the exponentiated one:

(x) = N i< j (x i -x j ) -→ N i< j 2 sinh x i -x j 2 . ( 36
)
In this case, insertion of an external source implies the expectation value of the Wilson loop operator, which is characterized by the Schur function [START_REF] Eynard | [END_REF].

To see this duality more explicitly, we now introduce another correlation function by multiplying the weight functions:

N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 = N ,M; p,q {a i } N i=1 , {b j } M j=1 ; {λ α } p α=1 , {μ β } q β=1 p α=1 e -1 W (λ α ) q β=1 e 1 W (μ β ) . ( 37
)
This can also be written in the determinantal form (31), but by replacing the auxiliary functions as follows:

P k (λ) → λ k e ∓ 1 W (λ) , Q k (a) → dx x k e ∓ 1 W (x)±xa , (38) R(a; b) → dxdy 1 x -y e -1 W (x)+ 1 W (y)+xa-yb , R(λ; μ) → 1 λ -μ e -1 W (λ)+ 1 W (μ) , (39) S L (a; μ) → dx 1 x -μ e -1 W (x)+ 1 W (μ)+xa , S R (λ; b) → dy 1 λ -y e -1 W (λ)+ 1 W (y)-yb . ( 40
)
The sign factor in (38) depends on whether M + p > N + q or M + p < N + q. In this case, it is easier to see the relation between auxiliary functions, which is summarized in Fig. 1.

PTEP 2014, 123A01

T. Kimura

Integrable equations

In the previous sections, we have shown that the partition functions of the Hermitian supermatrix models have a determinantal formula, and such a determinant structure plays an important role in the relation to the integrable systems [6]: it behaves as a τ -function for the corresponding integrable equation. In this section we study the connection between the supermatrix partition function and the integrable system.

Supermatrix model with external source

We now show that the supermatrix integral involving the external source ( 12) can be interpreted as the τ -function for the Toda lattice equations. In this paper we especially derive the integrable equation, obtained by taking all the external source parameters to the same value. We first split the external source parameters into the center of mass and relative parts:

a i = δa i + a, b j = δb j + b, i = 1, . . . , N j = 1, . . . , M . ( 41
)
With this separation, the determinant ( 7) is given by

N ,M (a; b) = N (δa) M (δb) N i=1 M j=1 (δa i -δb j + a -b) -1 . ( 42
)
Let us then consider the Taylor expansion around the center of mass for the functions ( 14) and ( 15):

Q k-1 (a i ) = ∞ l=1 (δa i ) l-1 (l -1)! Q k+l-2 (a), (43)
R(a i ; b j ) = ∞ l,m=1 (δa i ) l-1 (l -1)! (δb j ) m-1 (m -1)! R (l-1,m-1) (a; b), (44)
with

R (i, j) (a; b) = ∂ i ∂a i ∂ j ∂b j R(a; b). ( 45
)
By taking the limit δa i = δb j = 0 for ∀i, j, the determinant in ( 12) is written

det Q k-1 (a i ) R(a i ; b j ) = det 1≤i, j≤N (δa i ) j-1 ( j -1)! • det 1≤i, j≤M (δb i ) j-1 ( j -1)! • det Q i+k-2 (a) R (i-1, j-1) (a; b) , ( 46
) PTEP 2014, 123A01 T. Kimura
where the range of indices is given by i = 1, . . . , N , j = 1, . . . , M, and k = 1, . . . , N -M in the determinants of the LHS and the third one in the RHS. Therefore, after some cancellation, we obtain the following formula:

Z N ,M {a i = a} N i=1 , {b j = b} M j=1 = c N ,M ZN,M (a, b), (47) ZN,M (a, b) = det Q i+k-2 (a) R (i-1, j-1) (a; b) , ( 48
)
with the constant

c N ,M = (a -b) N M G(N + 1) -1 G(M + 1) -1 , (49)
where

G(n) = n-2 i=0 i! = n-1 i=1 (i)
is the Barnes G-function. We remark that this factor corresponds to the volume element for the U(N ) and U(M) groups. This shows that the supermatrix integral ( 4) is finally written as a kind of Wronskian in the equal parameter limit, δa i = δb j = 0.

In order to derive the integrable equation for the supermatrix integral, we now apply the Jacobi identity for a determinant

D • D i j k l = D i k • D j l -D i l • D cj k , ( 50
)
where D is a size-n determinant, and the size-n -1 minor determinant D i j is obtained by removing the ith row and the jth column from the matrix. Similarly, D i j k l is the size-n -2 determinant obtained by getting rid of the i, jth rows and the k, lth columns. Setting i = k = N and j = l = N -1, we have

ZN,M • ZN-2,M-2 = ZN-1,M-1 • ∂ a ∂ b ZN-1,M-1 -∂ a ZN-1,M-1 • ∂ b ZN-1,M-1 , (51)
and by putting i = k = N -M and j = l = N -M -1, we obtain

ZN,M • ZN-2,M = ZN-1,M • ∂ 2 a ZN-1,M -∂ a ZN-1,M 2 . ( 52
)
These yield the Toda lattice equations [START_REF] Hirota | The Direct Method in Soliton Theory: Cambridge Tracts in Mathematics[END_REF]:

• 2D Toda equation: ZN+1,M+1 • ZN-1,M-1 ZN,M 2 = ∂ 2 ∂a∂b log ZN,M . ( 53
)
• 1D Toda equation:

ZN+1,M • ZN-1,M ZN,M 2 = ∂ 2 ∂a 2 log ZN,M . ( 54
)
This means that the partition function is the τ -function for the Toda lattice system both in one and two dimensions simultaneously. This property can also be found in the ordinary Hermitian matrix model [5].

Characteristic polynomial with external source

We then show that the characteristic polynomial average with an external source (31) similarly obeys the integrable equations. In addition to the external source, we split the parameters for characteristic PTEP 2014, 123A01

T. Kimura polynomials:

λ α = δλ α + λ, μ β = δμ β + μ, α = 1, . . . , p β = 1, . . . , q . ( 55
)
As in the previous case, by considering the Taylor expansion around the center of mass, we obtain a Wronskian-type determinantal formula in the equal parameter limit δa i = δb j = δλ α = δμ β = 0 for ∀i, j, α, β, N ,M; p,q (a, b; λ, μ) = c N ,M; p,q ˜ N ,M; p,q (a, b; λ, μ), (56)

˜ N ,M; p,q (a, b; λ, μ) = det ⎛ ⎜ ⎝ Q i+k-2 (a) P (α-1) k-1 (λ) R (i-1, j-1) (a; b) S (α-1, j-1) L (λ; b) S (β-1,i-1) R (a; μ) R(α-1,β-1) (λ; μ) ⎞ ⎟ ⎠ , ( 57
)
with the constant

c N ,M; p,q = (a -b) N M (λ -μ) pq G(N + 1)G(M + 1)G( p + 1)G(q + 1) . ( 58
)
Here

P (α) k (λ) = d α dλ α P k (λ), S (α, j) L (λ; b) = ∂ α ∂λ α
∂ j ∂b j S L (λ; b), and so on. In this case, since the determinantal formula (57) consists of six blocks, we correspondingly obtain six equations by applying the identity (50):

˜ N ,M; p,q • ˜ N ,M; p-2,q-2 = ˜ N ,M; p-1,q-1 • ∂ λ ∂ μ ˜ N ,M; p-1,q-1 -∂ λ ˜ N ,M; p-1,q-1 • ∂ μ ˜ N ,M; p-1,q-1 , ( 59
) ˜ N ,M; p,q • ˜ N -2,M; p,q-2 = ˜ N -1,M; p,q-1 • ∂ a ∂ μ ˜ N -1,M; p,q-1 -∂ a ˜ N -1,M; p,q-1 • ∂ μ ˜ N -1,M; p,q-1 , ( 60
) ˜ N ,M; p,q • ˜ N ,M-2; p-2,q = ˜ N ,M-1; p-1,q • ∂ λ ∂ b ˜ N ,M-1; p-1,q -∂ λ ˜ N ,M-1; p-1,q • ∂ b ˜ N ,M-1; p-1,q , (61)
˜ N ,M; p,q • ˜ N -2,M-2; p,q

= ˜ N -1,M-1; p,q • ∂ a ∂ b ˜ N -1,M-1; p,q -∂ a ˜ N -1,M-1; p,q • ∂ b ˜ N -1,M-1; p,q , ( 62
) ˜ N ,M; p,q • ˜ N ,M; p-2,q = ( p -1) ˜ N ,M; p-1,q • ∂ a ∂ λ λ ˜ N ,M; p-1,q -∂ λ ˜ N ,M; p-1,q • λ∂ a ˜ N ,M; p-1,q , ( 63
) ˜ N ,M; p,q • ˜ N -2,M; p,q = p ˜ N -1,M; p,q • λ∂ 2 a ˜ N ,M; p-1,q -∂ a ˜ N ,M; p-1,q • λ∂ a ˜ N ,M; p-1,q . ( 64
)
We have used the following identity to derive some of these equations:

det ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ (x N ) (0) • • • (x N ) (M-1) . . . . . . (x N +M-2 ) (0) • • • (x N +M-2 ) (M-1) (x N +M ) (0) • • • (x N +M ) (M-1) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = M x det 1≤i, j≤M (x N +i-1 ) ( j-1) , (65)
PTEP 2014, 123A01

T. Kimura where we denote (x j ) (k) = d k dx k x j . From the first equation ( 59), e.g., we obtain the 2D Toda lattice equation

˜ N ,M; p+1,q-1 • ˜ N ,M; p-1,q-1 ˜ N ,M; p,q = ∂ ∂λ ∂ ∂μ log ˜ N ,M; p,q . ( 66
)
We can similarly obtain the 2D integrable equations from the others, except for (64), which leads to the 1D Toda lattice equation. Therefore, in this sense, the characteristic polynomial average with the supermatrix model plays the role of the τ -function for the Toda lattice system. Let us comment on higher-order integrable equations, which the matrix partition function would satisfy under generic parametrization of the source. In order to derive such equations, it is convenient to introduce the Miwa coordinate to the time variables,

t n (A) = 1 n Tr A -n . ( 67
)
In this sense, since there are two and four kinds of source parameters in (48) and (57), we can correspondingly apply two and four series of the time variables, respectively. Although these variables are decoupled in the equal parameter limit, and we obtain the well known Toda lattice equations, they might interact with each other when we consider more generic parametrization. If so, it would be interesting to study such a situation providing more involved integrable equations.

Discussion

In this paper we have studied the supermatrix model involving the arbitrary potential function W (x) with the external source term. We have derived the determinantal formula for the partition function and the characteristic polynomial expectation value with this supermatrix model. Based on this formula, we have exhibited the duality between the external source and the characteristic polynomial, and then pointed out that this duality is just interpreted as a Fourier transformation. We have also shown that the partition function and the characteristic polynomial average satisfy the Toda lattice equation both in one and two dimensions, especially in the equal parameter limit. This implies that we can obtain the corresponding τ -function as the supermatrix integral with an external source as well as the ordinary Hermitian matrix model. As pointed out in Sect. 3, we can obtain various situations from the supermatrix model by taking the corresponding limit. From this point of view, it is interesting to study not only simple reductions of the characteristic polynomial average (23), e.g., M = 0, q = 0, but also the analytic continuation to negative numbers. Actually, the U(N |M) supermatrix integral is, at least perturbatively, equivalent to the non-supersymmetric matrix model through the analytic continuation M → -M, but with the two-cut solution, which breaks the originally symmetry, U(N + M) → U(N ) × U(M). This kind of equivalence has recently been gathering a great deal of attention in the area of string/M-theory [22]. For example, it would be interesting to see how the duality relation is affected by such an analytic continuation. It is expected that we can obtain more various kinds of correlation functions, and all the reductions could preserve the determinantal structure [START_REF] Strahov | [END_REF]28], through the analytic continuation.

From the random matrix theoretical point of view, it would be interesting to extend the present result to other situations, e.g., the OSp(N |M) supermatrix model, chiral ensemble, complex matrix model, two-matrix model, and so on. To investigate the matrix model with an external source, the HCIZ formula plays an important role in integrating out the angular part of the matrix. Thus it seems better to start by deriving the corresponding formula, e.g., by considering the differential equation for the matrix integral.

Fig. 1 .

 1 Fig. 1. Fourier transform web: the relation between the auxiliary functions. FT x means Fourier transformation with respect to x variable. The duality (35) is seen as a double Fourier transformation with λ and μ variables.

2/12at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

3/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

5/12at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

7/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

8/12at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

9/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

10/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

11/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

12/12 at CEA on December 10, 2014 http://ptep.oxfordjournals.org/ Downloaded from

Acknowledgements

We would like to thank B. Eynard for the collaboration on related subjects. We are also grateful to M. Bergère and E. Brézin for correspondence. This work is supported in part by a Grant-in-Aid for JSPS Fellows (#25-4302).


Adam Gouge

Erwan Bocher

Nicolas Fortin

Gwendall Petit

A tool for understanding the influence of urban mobility plans (UMP) on spatial accessibility

Keywords: Urban mobility plans (UMP), spatial accessibility, shortest paths

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

Introduction

The Law on air and the rational use of energy (LAuRE), published on December 30, 1996, states that everyone has the right to breathe air that does not harm their health. In particular, it requires cities of more than 100,000 inhabitants to establish an Urban Mobility Plan (UMP1 ) in order to reduce the impact of public transportation on the environment. This law also requires the UMP to be reevaluated every five years, though it does not specify how.

In collaboration with Nantes Métropole, the Institut de Recherche en Sciences et Techniques de la Ville (IRSTV 2 ), the French Research Institute in Urban Sciences and Techniques, initiated a research project known as Eval-PDU (2008[START_REF] Bocher | OrbisGIS: Geographical information system designed by and for research[END_REF] in an effort to evaluate the environmental impacts of public policies on public transportation. The methodology includes numerous approaches, including numerical simulations evaluating air quality, noise pollution, energy consumption, and greenhouse gas emissions with a chain of physically-based models, based on alternative and comparative scenarios.

In this paper, we examine how to evaluate UMPs geographically. In effect, changing the kinds of transportation offered to the public (e.g., introducing a new tram line or installing a public bicycle system) or modifying the structure of the road network to favor public transportation (pedestrianizing a street, introducing bus lanes) influences how people choose to get around on a daily basis. Everyday activities such as going to work can quickly become a burden if the changes introduced by the UMP are not sufficiently evaluated and their impacts anticipated (the distance from point A to point B could change; access to public services could decrease).

To address these issues, we propose a methodology and an open-source tool which, starting from geographic data representing the public transportation network of a city, allows one to produce indicators representing accessibility to principle destinations, such as train stations, airports, universities, hospitals, gyms, etc. We further clarify these accessibility indicators by crossing them with demographical data, thus quantifying the impact of a UMP on a city's inhabitants by studying their spatial distribution relative to principal destinations.

Methodology

The very definition of a UMP, by modifying public transportation options available to citizens as well as traffic flow, implies a restructuring in the way roads are shared. This reorganization has direct repercussions on citizens' transportation choices around and, by extension, on their quality of life. When we realize how much time people spend in public transportation-an hour and a half per day on average for inhabitants of the Ile de France region; over two hours for 20% of them-the vital necessity of tools for quantifying and characterizing the impact of UMPs becomes clear.

A number of travel-related problems can be cast in the language of graph theory, which examines relationships among discrete entities. The entities are represented as a set V of vertices or nodes and the relationships as a set E of edges [START_REF] Newman | The structure and function of complex networks[END_REF]. The resulting graph is the ordered pair G = (V, E). In a travel context, the edges are often road segments and the nodes are intersections or points of interest.

A variety of techniques and indicators for studying the structure of the graph are well-known. Some examples include degree distributions, clustering, network correlations, random graph models, shortest path calculations, and closeness and betweenness centrality. [START_REF] Rodrigue | The geography of transport systems, Routledge[END_REF][START_REF] Newman | The structure and function of complex networks[END_REF]. In particular, Dijkstra's algorithm provides the solution to finding the most efficient travel route and the information that can be concluded from this calculation, such as identifying the closest facility or defining service areas based on travel time [START_REF] Rodrigue | The geography of transport systems, Routledge[END_REF]). Yet these results are not widely implemented in Geographical Information Systems (GIS) or other cartographic contexts.

In this paper, we examine the notion of accessibility, defined as "the measure of the capacity of a location to be reached by, or to reach different locations. Therefore, the capacity and the structure of transport infrastructure are key elements in the determination of accessibility" [START_REF] Rodrigue | The geography of transport systems, Routledge[END_REF]). In our case, we study accessibility via three different H2Network functions:

1. ST Graph: From a table of geographical data, we construct a mathematical graph by assigning numerical ids to edges and nodes, resulting in two new tables. An edge may be enriched with additional information that controls how shortest paths are calculated, such as its orientation or weight. For example, one can specify whether a road is one-way or bidirectional, as well as its length or estimated travel time. 2. ST Accessibility: We identify a subset of destination nodes (such as all hospitals). Fix any node v. Calculate the shortest distance from v to each of the destination nodes. Take the minimum of these distances, and designate the corresponding destination as the closest destination. Repeat this for every node v in the graph. 3. ST ShortestPathTree: Calculate the part of the graph which can be reached starting from a given start node and limited by a given time or distance.

We use the information obtained during these calculations to produce a service area map, for various services such as train stations or schools. This map includes isolines whose intervals are expressed in metric or temporal units. Finally, the isoline map is crossed with demographical data in order to quantify the number of inhabitants affected by changes in the road network introduced by a UMP. Using this information, we can estimate, for example, the percentage of inhabitants situated at a distance of less than 500 meters from a train station, or at a travel time of greater than 20 minutes from a school.

Implementation in H2Network

The methodology described above is implemented in the open-source Geographical Information System OrbisGIS [START_REF] Bocher | OrbisGIS: Geographical information system designed by and for research[END_REF], which uses the relational database management system (RDMS) known as H23 . In 2013, the OrbisGIS team made the choice to switch to H2 from the previous database they had invented, called GDMS [START_REF] Bocher | GDMS: an abstraction layer to enhance spatial data infrastructures usability[END_REF]. This choice brought greater modeling flexibility, improved performance and better interoperability with other databases. When the choice was made, H2 did not yet have a Geometry data type. The team collaborated with the H2 team to add basic support for Geometries, and eventually created their own open-source spatial extension of H2 called H2GIS4 .

H2GIS contains multiple spatial functions, including those specified by the OGC's Simple Features for SQL (SFSQL) 1.2.1 standards [START_REF] Herring | OpenGIS Implementation Specification for Geographic information -Simple feature access -part 2: SQL option[END_REF][START_REF] Herring | OpenGIS Implementation Specification for Geographic information -Simple feature access -Part 1: Common architecture[END_REF]. It includes additional spatial functions the team has developed that extend and complement the SFSQL functions. Examples include:

• ST CompactnessRatio, which computes the perimeter of a circle whose area is equal to a given polygon's area, and returns the ratio of this computed perimeter to the polygon's perimeter (useful for urban morphology), and • ST Delaunay and ST ConstrainedDelaunay, which calculate (constrained) Delaunay triangulations (useful in hydrology).

Both H2 and H2GIS are fully written in Java. They can be used in embedded mode in order to give spatial database capabilities to any Java application with a small mem-ory footprint. Moreover, H2GIS is based on a robust and well-known open source library known as the JTS Topology Suite5 , which provides an API for modeling and manipulating 2-dimensional linear geometry, as well as numerous geometric predicates and functions.

H2Network is the network analysis component of H2GIS. It provides functions for constructing a mathematical graph from geographical data representing a road network and performing shortest path calculations and other network analysis on the resulting graph. Some of its functions were mentioned in the previous section.

Case study

We take as study area the 24 municipalities of Nantes Métropole (566 500 inhabitants). The road network is modeled as 53,616 segments (polygonal chains); the data comes from the national topographical database BD TOPO R produced by the Institut National de l'Information Géographique et Forestière (IGN6 ), the French National Geographic Institute. The study area is divided into 229 IRIS zones, which are defined by the Institut national de la statistique et des études économiques (INSEE7 ), the French national statistics institute, to be statistical units comprising a population of 2,000 inhabitants. We represent them as polygons.

We study accessibility with respect to three scenarios:

1. Base case: The original road network as defined by the BD TOPO R data. 2. Busway: Scenario 1 with a new central bus line (the Busway) crossing Nantes from north to south and reducing light vehicles' access to the city center. 3. Bridge: Scenario 2 with a new bridge (the Eric Tabarly Bridge) which joins two major neighborhoods of Nantes, Ile de Nantes and Malakoff, otherwise naturally separated by the Loire river.

Overall, our results show a spatial redistribution of distances, which are notably more pronounced in Scenario 2. In effect, the accessibility calculations for light vehicles (cars) in this scenario show an increase in distances when approaching the center, as compared to Scenario 1. If we relate these distances to the inhabitants, our analysis shows that the number of inhabitants affected is relatively more important due to the fact that neighborhoods near the city center are more densely populated.

Invitation

We would like to invite you to attend and participate in our workshop entitled Road network analysis with H2Network: Applications of the spatial database H2GIS. We will demonstrate and apply the various analyses evoked in this paper.

http://www.installationsclassees.developpement-durable.gouv.fr/Urban-Travel-Plan-PDU.html, retrieved April

, 2014. 2 http://irstv.fr, retrieved April, 2014.

www.ogrs-community.org

http://www.h2database.com, retrieved April, 2014.

http://www.h2gis.org, retrieved April, 2014.OGRS 2014

http://sourceforge.net/projects/jts-topo-suite, retrieved April, 2014.

http://www.ign.fr/, retrieved April, 2014.

http://www.insee.fr, retrieved April, 2014.www.ogrs-community.org

http://belgrand-gebd.ifsttar.fr, retrieved April, 2014.OGRS 2014

Acknowledgments

The H2GIS team would like to thank the French GEBD-Belgrand 8 project for supporting this research.