Martin Golasowski
email: martin.golasowski@vsb.cz

Radek Tomis
email: radek.tomis@vsb.cz

Jan Martinovič
email: jan.martinovic@vsb.cz

Kateřina Slaninová
email: katerina.slaninova@vsb.cz

Lukáš Rapant
email: lukas.rapant@vsb.cz

Performance Evaluation of Probabilistic Time-Dependent Travel Time Computation

Keywords: probabilistic time-dependent route planning, Monte Carlo simulation, probabilistic speed profiles, Xeon Phi, HPC

Computational performance of route planning algorithms has become increasingly important in recent real navigation applications with many simultaneous route requests. Navigation applications should recommend routes as quickly as possible and preferably with some added value. This paper presents a performance evaluation of the main part of probabilistic time-dependent route planning algorithm. The main part of the algorithm computes the full probability distribution of travel time on routes with Monte Carlo simulation. Experiments show the performance of the algorithm and suggest real possibilities of use in modern navigation applications.

Introduction

Nowadays, there are many navigation services that help to find the best route to a requested destination. Drivers usually prefer the fastest route to the destination, which naturally depends on a departure time. Navigation services can find a variety of routes depending on many factors. Such factors can be map data, an algorithm used for the route planning, a possibility to use traffic data, etc. Recently, the traffic data are used extensively, because more precise results can be produced with their help. Many navigation services use traffic data to produce better and more driver friendly routes. However, such navigation services usually use only averages [START_REF] Rice | A simple and effective method for predicting travel times on freeways[END_REF] of historic travel times on roads for computation of the routes, even though this is not sufficient in many cases. Traffic events can significantly affect the actual travel time of a route [START_REF] Fan | Arriving on time[END_REF][START_REF] Hofleitner | Learning the dynamics of arterial traffic from probe data using a dynamic bayesian network[END_REF]. There can be recurrent traffic congestions at specific times which causes delays in traffic. Navigation services then should calculate with such events [START_REF] Miller-Hooks | Adaptive least-expected time paths in stochastic, time-varying transportation and data networks[END_REF]. The probability of the traffic congestion can be low, but the delay in the case the congestion happens can be very long. Therefore, it is advantageous to count with uncertain traffic events and their probabilities in the computation of the route [START_REF] Nikolova | Optimal route planning under uncertainty[END_REF][START_REF] Yang | Multi-cost optimal route planning under time-varying uncertainty[END_REF].

Our main motivation for this paper is to give drivers some added value for suggested routes. This added value is a complete probability distribution of travel time on a suggested route. The probability distribution of travel time on a route is a valuable information, because it indicates how likely it is for the driver to reach the destination according to travel time. We use our previously developed algorithm [START_REF] Tomis | Probabilistic time-dependent travel time computation using monte carlo simulation[END_REF] based on Monte Carlo simulation for the computation of probability distributions of travel time. The algorithm is briefly described at Section 2.

The added value in the form of the probability distribution of travel time brings increased computational complexity. As the first step, we suggest the best route like other navigation services, and as the second step we compute the probability distribution of travel time on this route, which requires additional computational time. The probability distribution of travel time on single route changes depending on a departure time and also on a day of the week. Traffic congestions are often connected to rush hours and some days of the week can have worse congestions than others, which is reflected in the probability distribution.

Many users might request a route at the same time and navigation services should be able to handle as many route requests at once as possible. There are two main ways to deal with a lot of route requests. The first one is to dynamically handle requests at the time they arrive and the second one is to perform the preprocessing of the probability distributions on all routes for all possible departure times and then handle requests from this precomputed data. However, it is almost impossible to cover all departure times in preprocessing, so some discretization of departure times have to be introduced. The computation time of such preprocessing is huge and depends on the number of departure times.

The main goal of this paper is to evaluate the limitations of our algorithm in a scenario with large amounts of simultaneous route requests in conjunction with a specific hardware. In other words, it would be beneficial to know how many route requests can be handled at once with different hardware. We should note that this paper evaluates only our previously developed algorithm [START_REF] Tomis | Probabilistic time-dependent travel time computation using monte carlo simulation[END_REF] for the computation of the probability distribution of travel time, because the regular route planning used for the selection of the route is an ubiquitous topic. The organisation of the paper is as follows. Section 2 describes the tested algorithm, Section 3 presents experiments and an architecture overview of a testing hardware. Finally, Section 4 concludes the paper with final remarks.

Algorithm Overview

This section presents a brief description of the algorithm [START_REF] Tomis | Probabilistic time-dependent travel time computation using monte carlo simulation[END_REF] used for the computation of the probability distribution of travel time. The input for the algorithm is some selected route and a departure time. The route is composed as a line of road segments. A Monte Carlo simulation is used for the computation of probability distributions. The simulation randomly selects probabilistic speed profiles on road segments and computes travel time at the end of the route. Probabilistic speed profiles are time-dependent sets of speeds and their probabilities assigned to road segments. Many Monte Carlo simulation iterations are needed to obtain enough travel times for the construction of the probability distribution of travel time. The number of simulation iterations greatly affects the precision of the result.

Probabilistic speed profiles present a way to state that various levels of traffic can occur with some probability. Different profile types can exist for variety of situations, e.g. free flow, congestions, traffic jams. All individual profiles are determined by belonging road segment and the time interval of its validity. Probabilistic speed profiles can be created in many ways. One way is to use Level-of-Service (LoS), which is the measured speed divided by the free flow speed. We use this method and describe it in detail at Section 3.

Algorithm 1 shows a pseudo code of Monte Carlo simulation. The simulation is repeated for a specified number of samples, and then the probability distribution is created from computed travel times. Time t is an actual travel time at the end of a road segment. Notice that the computation of t for a particular road segment depends on t from its preceding road segment. This makes single simulation computationally dependent and hard to parallelize. However, individual simulations are very fast and strictly independent of each other and can be run in parallel. for s in segments do 5: r = random number between 0 and 1 6: t = ComputeTravelTime(s, t, r) speed profile is selected by r 7:

end for 8:

add travel time t to result 9:

end for 10: end procedure

Experiment

Our algorithm is based on the well-known Monte Carlo method which is often regarded as embarrassingly parallel algorithm. The computation involves running a large number of independent instances of the algorithm with different data which makes it perfectly suitable for massively parallel computation platform such as Intel Xeon Phi. However, various problems can arise since sequential part of our implementation is memory/branch bound and is unable to fully exploit the power of vector processing unit available on the coprocessor. We expect good scalability on the coprocessor as well as on the host CPU even with this potential handicap. In this section we present results measured on the coprocessor and optimizations applied to our code to gain better speedup.

We do not try to establish a benchmark of the coprocessor since our code is not a representative sample of all Monte Carlo implementations. However, every optimization applied can improve run time on the coprocessor as well as on the host CPU. Therefore, at the end of this section we present a comparison of run times on both platforms.

Intel Xeon Phi Coprocessor Architecture Overview

The Intel Xeon Phi is a SIMD oriented coprocessor based on the many integrated core (MIC) architecture. The first version known as Knights Ferry was introduced in 2010 as a proof-of-concept targeted at high performance computing segment to counter-balance increased demand for GPGPU based architectures. The Phi has 61 x86-based cores, each core capable of running 4 threads at a time [START_REF] Rahman | Intel R Xeon Phi TM Coprocessor Architecture and Tools: The Guide for Application Developers[END_REF]. The cores are based on heavily modified P54c architecture extended by new vector processing unit (VPU) and larger cache memory. The cores communicate with the memory via ring interconnect bus based on the uncore topology introduced by Intel with the Sandy Bridge microarchitecture. It has 8 distributed memory controllers communicating with graphics double data rate (GDDR5) RAM memory. The connection with host processor is realized via 16x PCI express bus.

The main power of the Xeon Phi is in the new VPU which supports instructions from the AVX-512 instruction set. Using these instructions, each core of the Xeon Phi is able to perform arithmetical operations on either 16 single precision or 8 double precision floating point numbers at once. A theoretical speed up of 16x is achievable by correct usage of the Xeon Phi VPU. However, correct usage of the VPU depends on various factors, such as data access pattern, memory alignment, and size of the data set (ideally, the entire working data set should fit in the on-board memory, since PCIe bus introduces a significant bottleneck).

Thanks to its roots in x86 architecture, the programming model for the Phi is more or less the same as for the regular CPUs. As the card itself runs a common Linux-based operating system, standard threading models such as OpenMP, Intel TBB, or pthreads can be used, as well as MPI. Intel supplies its own SDK with optimized version of the math kernel library (MKL) and C/C++ and Fortran compilers. However, to achieve optimal exploitation of the Xeon Phi computational power, a substantial amount of attention should be paid to effective partitioning of the workload and to the overall differences of the architecture by profiling and optimization of the code.

Although cores of the Xeon Phi have superscalar properties, the cores itself are relatively slow (1.2 GHz < 2.5 GHz) compared to the Haswell CPUs and lack several features such as level 3 cache or out of order execution. Due to this deficiency, the Phi cores may perform worse in sequential or branch-heavy code compared to the host CPU (such as Haswell). On the other hand, well vectorized code can run significantly faster on the Phi than on the host CPU, thanks to the VPU and significantly larger (12 < 61) number of available cores.

Traffic Data

The data for our experiment come from motorways and 'A' roads managed by the Highways Agency, known as the Strategic Road Network, in England 12 . This data contain time series of average journey time, speed and traffic flow information for segments of various length (from few hundred meters to several kilometers) with 15-minute aggregation periods (i.e. four values for each hour). Traffic speed in these data sets is calculated using a combination of stationary sensors (Automatic Number Plate Recognition, Induction loops) and floating car data sources.

In our experiment, we used data from January to December 2014. Route from Sunderland to Thetford was chosen, because it is one of the longest direct routes contained in the data set. This data were used to calculate probabilistic speed profiles for the segments along the route and to calculate empirical free flow speed for these segments. These free flow speeds were calculated as average of top 10% aggregated speed measurements on the chosen segment (as these speeds are most probably not influenced by the traffic density).

Profiles used for experiments were created for four LoS intervals, specifically 100; 75), 75; 50), 50; 25) and 25; 0 . An example of such speed profiles for the entire route from Sunderland to Thetford, Monday 8:00 is presented in Figure 1.

The interpretation of what these intervals mean can be as follows. The first interval 100; 75) represents free traffic, the second interval 75; 50) corresponds to light traffic, the third interval 50; 25) means heavy traffic, and the last one 25; 0 is for traffic jams.

Experimental Simulations

The simulation was performed with the traffic data described in Subsection 3.2 on the Salomon supercomputer operated by IT4Innovations National Supercomputing Centre, the Czech Republic. Nodes of the Salomon cluster contain a of the Xeon Phi coprocessor based on the Knights Corner microarchitecture (7210P). This model has 16 GB of on board RAM memory, 61 compute cores, and there are 2 cards per compute node of the cluster. In total, there are 432 nodes containing 2 Xeon Phi coprocessors. Each compute node has two 12 core Haswell CPUs and 128 GB RAM.

An example of simulation output is presented in Figures 3a and3b. The single profile simulation (a) uses a selected profile as a starting point, then runs the simulation for the entire path. In this case, the starting point was determined by week day (Tuesday) and time (8:00 AM). The all profiles simulation (b) runs the algorithm for all available speed profiles determined by its parameters. (i.e. days of the week and 15 minute intervals of a single day). Both the simulations have been run with 1000 Monte Carlo samples. The first simulation provides information about the most probable travel time along the given path and for a selected time of departure. Based on results of our simulation shown in Figure 3a, it can be said that if you drive from Sunderland to Thetford on Tuesday around 8:00 AM, you will most probably reach your destination in about 6 hours and 30 minutes (23,200s). By looking at results of the simulation for all profiles in Figure 3b, it is evident that the estimated travel time can vary through the week, however the journey will most likely take no more than 5 hours and 30 minutes (19,800s).

Performance on the Intel Xeon Phi Coprocessor

Various other factors involving configuration of the run time environment can affect performance of the coprocessor. One of the factors in our case was thread allocation strategy. Our code is parallelized using the OpenMP standard. Intel provides its own OpenMP implementation as part of the Xeon Phi Manycore Software Stack (MPSS) toolchain. The library allows to choose between various thread allocation strategies by setting KMP AFFINITY environment variable .

There are three main strategies. The compact strategy allocates whole cores in blocks. For example, running 40 threads on the Phi would use only 10 cores since single core runs 4 hardware threads. The scatter strategy assigns threads to cores evenly in round-robin fashion across all cores. The explicit strategy is the last one and allows explicit assignment of a given software thread to a particular core.

Measurements of the code scalability with the compact and scatter strategy are presented in Figure 4. The scatter strategy is more suitable in our case since our code is memory/branch bound and running more than 2 threads per core creates unwanted congestion during memory reads. Another factor involved is the cache hierarchy of the coprocessor. The Xeon Phi has 32kB + 32kB (instruction, data) L1 cache and 512kB L2 data cache per core. The graph in Figure 5 shows differences in performance of the code on cold and hot caches. The cold cache times were measured from single execution of the code on the coprocessor for different number of threads. The hot cache times were measured as average of 10 successive runs of the code.

The graphs show that proper use of the cache memory can provide significant speed up of our code. Unfortunately, the input data for single run will most likely differ each time, since most of the user requests will be computed for different road segments and start time. The real world usage of the coprocessor is reflected by the cold cache measurement which means we are able to effectively use only half of the resources available for efficient usage of the coprocessor.

Scalability of the code on the host CPU is presented in Figure 6. The code scales very well while execution time is approx. 2 times shorter when running on single 12 core Haswell CPU that on the coprocessor. The code running on the CPU benefits from L3 data cache as well as from out-of-order execution and more sophisticated memory prefetching. Measured times of execution on the coprocessor are presented in the Table 2 and measured times on the host CPUs in Table 3. All measurements were performed with 1,000 samples run, threads were distributed evenly across all available cores of the target platform (except for the compact thread affinity strategy).

Code Optimization

In parallel Monte Carlo based simulations, the optimization effort usually focuses on finding a suitable method for random number generation and on optimizing data access pattern since only single synchronization occurs at the end of all runs and no data dependencies exist between the threads. The core of our algorithm relies on the selection of a speed profile given a canonical probability p ∈ (0; 1). Effective indexing of the speed profiles for a given path segment can significantly improve overall run time of the algorithm. The speed profile for single time interval (identified by a day of the week and a time) is represented by pairs of velocity and probability values (see Table 3.2). The velocity values are distributed in an array of fixed length according to their probability. The size of the array determines the resolution of the probability values which can be represented by a linear index. An example of the mapping between speed profiles and the storage array is presented in Table 4. The first table shows a single probabilistic speed profile with four LoS. The second table shows the final distribution and encoding of the velocities to the linear array indexes according to their probability. All the speed profiles for single road segment are represented by matrix of the velocity arrays stored in single linear array to improve usage of the cache memory. This type of the representation lowered the amount of necessary processor instructions used for the selection of the appropriate speed from the speed profile. The early version of our code computed the index of the profile in a given day from time specified in UNIX timestamp format. The computation of the profile index involved several system function calls (namely localtime) which posed a significant bottleneck in the multithreaded environment on the coprocessor. The next version was based on simple data types and did not use any date-time related functions from the C standard library.

However, profiling in Intel VTune 4 suggested that another bottleneck appeared as we used pseudo-random number generation routines (introduced in C++11) available in C++ Standard Library. The slowdown was apparent especially when using large number of threads such as on Xeon Phi accelerator, since the library implementation was mainly sequential and thus unsuitable for usage on the Xeon Phi due to lack of any parallelism and vectorization.

Intel provides optimized version of the math kernel library (MKL) 5 for the Xeon Phi in its SDK 6 and we have decided to use their Mersenne-Twister implementation. This implementation uses vector instructions and allows us to generate the random numbers in batches which leads to more efficient use of the coprocessor. In the current version of our code, the random number generation takes only approx. 20% of the computation time according to the Intel VTune profiler in comparison with 90% in the previous version of the code.

Conclusion

The paper evaluated the computational performance of our algorithm for the computation of probabilistic time-dependent travel time on the Intel Xeon Phi coprocessor and Haswell CPU. Various code optimizations described in Section 3.5 have been done with appropriate speed up. We expected good scalability on the Xeon Phi coprocessor due to large parallelization potential of the algorithm even though it was not taking any advantage of the vector processing unit, and seemed generally unsuitable for the coprocessor. However, we have been able to use the VPU by using Intel MKL Vector Statistics Library. Finally, the good scalability was achieved on the coprocessor as well as on the host CPU. Measured execution times and scalability graphs are presented in Section 3.3. By performing experiments, we estimated that it is the best to run our code on 122 threads evenly distributed across all 61 cores of the coprocessor.

The experiments met our expectations, however a number of specific optimizations had to be applied to obtain these results. Each optimization that led to the improvement for the coprocessor also improved the execution time on the CPU, often significantly. Finally, the execution time was at best 2 times slower on the coprocessor than on the single host CPU. In general, we have been able to lower the overall execution time quite significantly. The early version of the code ran for approx. 5 minutes on single coprocessor while optimized version of our code runs for approx. 400 ms for the same input data.

It follows that the algorithm can easily handle many simultaneous route requests and the number of requests should be large enough to utilize the full power of the available hardware. Using the current implementation we can offload half of the work to the coprocessor and efficiently utilize the entire power of the HPC infrastructure.

In the future work, we would like to focus on further optimizations related to the cache memory subsystem of the coprocessor and to finding a vectorizable version of the algorithm. Another possible area of the research emerged during the optimization process. Some of the optimization techniques can be generalized for this type of Monte Carlo based algorithms by using some type of domain specific language, for example DSL and autotuning techniques being developed within the ANTAREX7 project [START_REF] Silvano | The antarex approach to autotuning and adaptivity for energy efficient hpc systems[END_REF]. This way, the specific optimizations can be applied by users without the need to perform such extensive studying of the target hardware microarchitecture.

Algorithm 1

 1 Monte Carlo Simulation for the Computation of the Probability Distribution of Travel Time for Given Array of Road Segments, Departure Time and Number of Samples

Fig. 1 .

 1 Fig. 1. Speed Profiles for the Selected Route, Monday 8:00 AM

Figure 2

 2 Figure 2 shows an example of speed profiles of the segments on the whole selected route for the interval 50-75% of LoS. The color indicates the probability that this interval of LoS occurrs. The speed profiles were created on basis of data from January to December 2014 considering days of the week. Table 3.2 shows example of probabilistic speed profiles for single road segment and for one hour. Each LoS is represented by two values, speed s in kilometers per hour and probability p ∈ 0; 1 .

Fig. 2 .

 2 Fig. 2. Heatmap of Speed Profiles for 50 -75% of LoS

Fig. 3 .

 3 Fig. 3. Example of the Simulation Output

Fig. 4 .

 4 Fig. 4. Scalability of Simulation for All Available Profiles on the Xeon Phi Coprocessor

Fig. 5 .Fig. 6 .

 56 Fig. 5. Scalability of Simulation for All Available Profiles on the Xeon Phi Coprocessor

Table 1 .

 1 Example of Probabilistic Speed Profiles

				Level of Service
	Day of Week Hour	100%	75%	50%	25%
			s [km/h] p [%] s [km/h] p [%] s [km/h] p [%] s [km/h] p [%]
		8:00 83.45 0.21 61.33 0.058 39.46 0.29 14.39 0.44
	Monday	8:15 83.40 0.19 56.25 0.086 36.18 0.21 15.44 0.52
		8:30 85.71 0.13 66.33 0.057 35.29 0.32 14.12 0.48
		8:45 84.00 0.17 63.08 0.230 40.23 0.25 17.05 0.34

Table 2 .

 2 Measured Run Times on the Coprocessor

	Threads	Times [ms] Scatter Compact Cold Hot
	1	40,111	40,207 40,162 40,118
	61	637	995	819	705
	122	410	571	731	514
	244	332	331	955	553

Table 3 .

 3 Measured Run Times on the Host CPU

	Threads	Times [ms]
	1	1,973
	12 (single CPU)	166
	24 (two CPUs)	85

Table 4 .

 4 Example of Probabilistic Speed Profile Encoding

		Speed	83.45 61.33 39.46 14.39
		Probability 0.21 0.058 0.29 0.442
	Index 0	1	2	3	4	5	6	7	8	9
	Speed 83.45 83.45 61.33 39.46 39.46 39.46 14.39 14.39 14.39 14.39

https://data.gov.uk/dataset/dft-eng-srn-routes-journey-times

http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/

https://software.intel.com/en-us/articles/how-to-analyze-xeon-phi-coprocessorapplications-using-intel-vtune-amplifier-xe-2015

https://software.intel.com/en-us/mkl-reference-manual-for-c

https://software.intel.com/en-us/articles/intel-manycore-platform-software-stackmpss

ANTAREX (AutoTuning and Adaptivity appRoach for Energy efficient eXascale HPC systems), Project Reference: 671623, Call: H2020-FETHPC-2014http://www.antarex-project.eu

Acknowledgment

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project 'IT4Innovations excellence in science -LQ1602', from the Large Infrastructures for Research, Experimental Development and Innovations project 'IT4Innovations National Supercomputing Center LM2015070', and co-financed by the internal grant agency of V ŠB -Technical University of Ostrava, Czech Republic, under the project no. SP2016/179 'HPC Usage for Transport Optimisation based on Dynamic Routing II'.


Lukashevich Maryna
email: lukashevich@bsuir.by

Starovoitov Valery
email: valerystar@mail.ru

An Approach to Cell Nuclei Counting in Histological Image Analysis

Keywords: histological analysis, machine analysis, cell nuclei counting

The paper describes a technique for automated cell nuclei counting. In this study, the primary goal is to provide simple and effective automated scheme of cell nuclei counting. The experiments on public data set of histology images have demonstrated acceptable level of calculation results.

Introduction and Motivation

The most important and rapidly developing areas in the field of visualization and control of microscopic objects is development of automated computer vision systems. Computer analysis of microscopy cell images has many real-life applications in a wide range of areas including diagnosis of a disease, morphological cell analysis and statistics [START_REF] Irshad | Methods For Nuclei detection, Segmentation, and Classification in Digital Histopathology: A Review Current Status and Future Potential[END_REF]. It plays an important role in biomedical research and bioinformatics. Microscopic analysis approach is particularly important in solving the following problems: malignant transformation of normal cells and cancer cell detection, morphological changes in cells, dynamic changes in the cells during therapeutic procedures [START_REF] Chen | Recent Advances in Morphological Cell Image Analysis[END_REF].

For that reason, automatic counting of cell nuclei is a key block in systems for microscopic analysis of cell images. Since 1970s automated methods of histological analysis have been developed [START_REF] Irshad | Methods For Nuclei detection, Segmentation, and Classification in Digital Histopathology: A Review Current Status and Future Potential[END_REF][START_REF] Chen | Recent Advances in Morphological Cell Image Analysis[END_REF]. Detection of cellular structures and cell counting are common tasks for many investigations. Some researchers have applied thresholding for cell counting in the investigated images. It is a simple way for separating objects of interest from the background. Another approach is to segment images using different edge detection (LoG filter, Laplacian filtering, etc.). However, detection and segmentation of cell nuclei is a challenging task, since the cells have a complex and a nonuniform structure. Another feature of histological image is non-uniform illumination, which also applies limitations on the use of standard approaches

In recent years many methods have been proposed for cell nuclei segmentation, separation and classification in histological analysis. More sophisticated approaches consist of several image processing stages. These methods are based on traditional image processing algorithms like adaptive contour model, watershed, morphological operations, k-means, Support vector machine, etc [START_REF] Jung | Impact of the Accuracy of Automatic Segmentation of Cell Nuclei Clusters on Classification of Cell Nuclei Clusters on Classification of Thyroid Follicular Lesions[END_REF][START_REF] Saharma | A Multi-resolution Approach for Combining Visual Information using Nuclei Segmentation and Classification in Histopathological Images[END_REF][START_REF] Alilou | Segmentation of cell nuclei in heterogeneous microscopy images: A reshapable temlates approach[END_REF][START_REF] Kowal | Nuclei Segmentation for Computer-Aded Diagmosis of Breast Cancer[END_REF][START_REF] Wienert | Detection and segmentation of cell nuclei in virtual microscopy images: a minimun-model approach[END_REF][START_REF] Zang | White Blood Cell Segmentation by Color-Space-Based K-Means Clustering[END_REF][START_REF] Song | Region-based progressive localization of cell nuclei in microscopic images with data adaptive modeling[END_REF][START_REF] Lurs | Murphy Nuclear segmentation in microscope cell images: a hand-segmented dataset and comparison of algorithms[END_REF][START_REF] Signolle | Wavelet-Based Multiscale Texture Segmentation in Application to Stromal Compartment Characterization on Virtual Slides[END_REF][START_REF] Lezoray | Marinette Revenu: Segmentation of cytological image using color and mathematical morphology[END_REF].

Recent works suggest combining different approaches for increasing performance. Several excellent reviews about methods for nuclei detection, segmentation and classification can be found in [START_REF] Irshad | Methods For Nuclei detection, Segmentation, and Classification in Digital Histopathology: A Review Current Status and Future Potential[END_REF][START_REF] Chen | Recent Advances in Morphological Cell Image Analysis[END_REF].

The main methods for cell nuclei segmentation, separation and classification are presented in Table 1. In our research we study the problem of automatically counting of cell nuclei in histological images of their size and shape.

Depending on the quality of prepared medication and used optical equipment, microscopy image can be both grayscale and color. The large number of different objects in histological images is a serious disadvantage for analysis of this class of images. Unfortunately, cell nuclei analysis is challenging task because of the complexity of the natural objects (noisy images, background heterogeneity, variations in object characteristics). Most of the cell nuclei are part of histological structures with complex and irregular visual features [START_REF] Loukas | An Image Analysis-Based Approach for Automated Counting of Cancer Cell Nuclei[END_REF]. Samples of microscopy images and their intensity histograms are presented in Fig. 1. Shape of histograms shows that the problem of automatic nuclei detection and segmentation in histological images is difficult by thresholding.

Methods and Materials

The investigated images are digital microscopy images stained with H&E (Hematoxylin & Eosin) from a publicly available data set [START_REF] Al-Kofahi | Improved automatic detection and segmentation of cell nuclei in histopathology images[END_REF][START_REF] Kofahi | Cell-basedquantification of molecular biomarkers in histopathology specimens[END_REF]. Authors [START_REF] Al-Kofahi | Improved automatic detection and segmentation of cell nuclei in histopathology images[END_REF][START_REF] Kofahi | Cell-basedquantification of molecular biomarkers in histopathology specimens[END_REF] used whole slide and tissue microarrays derived from tissue samples available through routine diagnostic. The histological slides were using the Zeiss Mirax Scan side scanner. The slide scanner was equipped with Zeiss Plan-Apochromat 20x (numerical aperture = 0.8) and an AVT Marlin F-146C Firewire 1/2 CCD camera with 4.65µm × 4.65µm pixel size. Combine with the 20x objective and 1x Cmount adapter the resulting image resolution is 0.23µm × 0.23µm. All slides were scanned at 20x and performed at full resolution. The resulting images were converted to the Virtual Slide Format with actual image data encapsulated and saved as JPEG image files with 85% JPEG quality. The image sizes vary from 33280 × 29184 pixels to 70080 × 159000 pixels, but the analyzed field of the images had a size of 600 × 600 pixels.

The final dataset contains 7931 cells from 36 color tissue images. The images of different organs (breast, kidney, gastric mucosa, connective tissue, small intestine, etc.) are found in data set. The nuclei center coordinates are used as ground truth data for proposed approach validation. All images were labeled by three pathologists, Fig. 2.

The Main Idea

The proposed approach for cell nuclei counting in histological image analysis consists of the following stages: 1. RGB color space regularization; 2. image intensity adjusting; 3. guided filtering; 4. thresholding; 5. morphological processing.

The preprocessing stage includes several steps. Original images contain a limited number of primary colors such as pink, purple and shades of blue. It is proposed to reduce the number of colors in the original image palette using Minimum Variance Quantization. Minimum variance quantization cuts the RGB color cube into smaller boxes (not necessarily cubical shape) of different sizes, depending on how the colors are distributed in the image. If the input image actually uses fewer colors than the number specified, the output colormap is also smaller, Fig. 3. For the tested data set this procedure is almost visible to the human eye because of the limited number of colors in the difference between original images and same images after regularization is hardly visible (Fig. 3), but it is crucial for further image analysis stages.

After that we try to increases the contrast of the images by a contrast enhancement procedure.

The next step is using guided filter. The filter performs edge-preserving smoothing of images, using content of a second image, called a guidance image, to influence the filtering. The guidance image can be the original image itself, a another version of the original image, or a completely different image. Guided image filtering is a neighborhood operation, like other filtering operations, but takes into account pixel statistics in the corresponding spatial neighborhood of the guidance image when calculating the value of the output pixel [START_REF] Xu | Image Smoothing via L0 Gradient Minimization[END_REF].

Then we apply binarization by Sauvola local image thresholding [START_REF] Sauvola | Adaptive document image binarization[END_REF][START_REF] Shafait | Efficient implementation of local adaptive thresholding techniques using integral images[END_REF][START_REF] Stathis | An Evaluation Technique for Binarization Algorithms[END_REF] for separation nuclei and background, see Fig. 4.

We clean the obtained binary images by opening and closing with a small structure element. After that we count automatically the number of black blobs corresponding to the cell nuclei. Opening and closing to the binary images we apply analysis help us to clean the binary images. They are both derived from the fundamental operations of erosion and dilation. Like those operators they are normally applied to binary images, although there are also gray level versions. The basic effect of an opening is somewhat like erosion in that it tends to remove some of the foreground (bright) pixels from the edges of regions of foreground pixels.

Nuclei Counting And Method Validation

We have tested implementations (in Matlab) of the proposed approach in 36 images from the data set described above.

The typical parameter values used in our experiments are shown in Table 2. Example of analysis results on H&E images is shown in Fig. 5. The performance of calculation was calculated using TP (true positive), FN (false negative) and FP (false positive) events. The performance parameters precision and recall calculating using (1)-( 2).

P recision =

T P T P + F P (1)

Recall = T P T P + F N (2)
A that combines precision and recall is the harmonic mean of precision and recall, the traditional F-measure or balanced F-score (3):

F = 2 • P recision • Recall P recision + Recall (3)
We compared the results with recent advances in this area [START_REF] Wienert | Detection and segmentation of cell nuclei in virtual microscopy images: a minimun-model approach[END_REF][START_REF] Song | Region-based progressive localization of cell nuclei in microscopic images with data adaptive modeling[END_REF][START_REF] Lurs | Murphy Nuclear segmentation in microscope cell images: a hand-segmented dataset and comparison of algorithms[END_REF], Table 3.

Discusion

Table 2 shows the best parameters for the tested image database. Developed approach can be adapted to different nuclei counting task. Its main advantage is automatically approach for cell nuclei counting. In future, it is necessary to automatization the process of the initial parametrs settings depending on images. Cell nuclei with small area were lost because we used the guided filter and mathematical morphology. But it allows us to get more accurate results for other cells. This is due to the fact that the improved thresholding results (disappear small areas that are not true cell nuclei). Future work will be focused on improving the segmentation stage and enlargement of the tested image data set.

In contrast to the known algorithms, the proposed approach does not segment all of the nuclei in images. The aim of the experiments was not to improve metrics Precision and Recall. We just calculated them with results of the similar research. In our method the value of F-measure is equal to 0.8743. It is bigger than F-measure of Al-Kofahi, but smaller than F-measure of Wienert results. There are a lot of small nucleus in the images, which have been marked by the experts. These nucleus difficult to correctly identify and count. In this way, we obtain decrease Precession accuracy. More subtle settings allow detecting the missing small nucleus, but this requires coordination with experts about nuclei features.

Conclusions

An automated cell nuclei detection and counting method based on the guided filter and morphological operations is proposed. Different to the existing algorithms, the proposed approach does not segment nuclei in images. It only detects the cell nuclei in a histological image and estimate the number of the detected nuclei. It avoids complicated algorithmic computations and provides good accuracy in the cell nuclei counting. The software prototype developed in our study may be considered as an automatic tool for a cell nuclei analysis.

Fig. 1 .

 1 Fig. 1. H&E image samples (a), (c), (e) and their intensity histogram (b), (d), (f) respectively.

Fig. 2 .

 2 Fig. 2. An example of H&E image with labeled nuclei centers.

Fig. 3 .

 3 Fig. 3. (a) An input image, (b) the image after requantization.

Fig. 4 .

 4 Fig. 4. Improvement in binarization by guided filtering: (a) input grayscale image, (b) thresholding without guided filtering; (c) input grayscale image after guided filtering; (d) thresholding after guided filtering.

Fig. 5 .

 5 Fig. 5. (a) An example of the original image and (b) its binary variant with the centers of cell nuclei labeled by experts.

Table 1 .

 1 Short Summary of State-of-the-Art Cell Nuclei Segmentation, Separation and Classifiction

	Solvable task	Methods
		adaptive contour model,
		adaptive thresholding,
		morphological operations,
	Segmentation	watershed, k-means,
		H-maxima transform,
		region growing,
		Gaussian mixture model
		distance transform,
	Separation	watershed transform, concavity detection,
		edge path selection
		texture,
		morphology,
		texture and morphology,
	Classification	texture/morphology/intensity with Support vector machine,
		texture/morphology/topology with Bayesian,
		intensity with k-means clustering,
		intensity&texture with Adaboost

Table 2 .

 2 Typical Parametrs Values

	Processing stage	Parameters	Different values
	Minimum Variance Quantization	amount of color space	4
	guided filtering	mask size balance between data matching and smoothing	33 0.5
	Sauvola local image	local thresholding with M-by-N neighbourhood	45×45
	thresholding	threshold	0.35
	morphology	radius of structuring elements shape of structuring elements	1 disk

Table 3 .

 3 Detection Results of the Proposed Approach, Wienert's, Al-Kofahis Methods in Precision and Recall Values

		Proposed	Wienert	Al-Kofahi
	Precision	0.833±0.01	0.908±0.04	0.707±0.13
	Recall	0.920±0.01	0.859±0.04	0.916±0.04
	F-measure	0.8743	0.8828	0.7980


Wojciech Bożejko
email: wojciech.bozejko@pwr.edu.pl

Łukasz Kacprzak
email: lukasz.kacprzak@pwr.edu.pl

Piotr Nadybski
email: nadybskip@pwsz.legnica.edu.pl

Mieczyslaw Wodecki
email: mieczyslaw.wodecki@uwr.edu.pl

Mieczysław Wodecki

Multi-machine Scheduling with Setup Times

de documents scientifiques de 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

In the vast majority of the considered in the literature scheduling problems, it is assumed that a task is done on a single machine. However, in many real production processes, in particular -modern computational systems -a task execution requires the use of more than one machine (CPU), Bożejko et al. [START_REF] Bożejko | Solving timetabling problems on GPU[END_REF]. Then, we can talk about systems with parallel machines and multi-machine (multiprocessor) tasks. We can meet such an approach also in cloud exploration [START_REF] Smutnicki | A new approach for multi-criteria scheduling[END_REF]. In literature, such problems have been considered for a long time, for instance in the work of Drozdowski [START_REF] Drozdowski | Scheduling multiprocessor tasks an overview[END_REF], [START_REF] Drozdowski | Scheduling for parallel Processing[END_REF] and in the monograph [START_REF] Drozdowski | Select problems of scheduling tasks in multiprocessor computer system[END_REF], Fleitelson [START_REF] Feitelson | Theory and practice in parallel job scheduling[END_REF], as well as in doctoral thesis of Kramer [START_REF] Kramer | Scheduling multiprocessor tasks on dedicated processors[END_REF] where they mainly concern computer systems.

The first applications of multi-machine models of tasks scheduling refer not only to chemical industry (Bozoki and Richard [START_REF] Bozoki | A branch-and-bound algorithm for thecontinous-process job-shop scheduling problem[END_REF]) but also to projects scheduling (Vizing [18]). The natural reference for this type of tasks are, examined from decades, multiprocessor computer systems (Błażewicz et al. [START_REF] Błażewicz | Scheduling multiprocessor tasks to minimize schedule length[END_REF]). In the vast majority of them the completion date of the execution of all tasks (C max ) was maximized. Today, many computer centers offer performance of calculations and the cost of the service depends not only on the duration of the computations but also on the number of required processors. Therefore, there appears a need for the construction and testing of new, corresponding to the reality models and analysis of criterion functions taking into account the overall costs. We can encouter similar issues in construction projects planning, where the simultaneous use of multiple resources is the basis for the organization of work.

The basic problem of multi-machine scheduling, P |size j |C max (where P -is a symbol of multimachine, whereas size j -required number of machines) was studied mainly in reference to the approximate algorithms and the worst case scenario(Lin [START_REF] Lin | Scheduling algorithm for nonpreemptive multiprocessor tasks[END_REF], Lloyd [START_REF] Lloyd | Concurrent task system[END_REF]). It has been proven that already in the case of two machines, the problem P 2|size j |C max is NP-hard. Błądek et al. [START_REF] Błądek | On contiguous and noncontiguous parallel task scheduling[END_REF] consider the problem with C max criterion assuming further that to certain tasks there must be assigned a fixed number of neighboring machines. In addition, it is assumed that the time of task scheduling undergoes learning or aging process, depending on the start of its execution.

In this paper we consider the problem of multi-machine scheduling, in which, between the tasks performed in sequence, there are additionally setups of machines introduced. In the case of construction processes, a setup is the time required for the movement of machines, equipment and employees. However, in computer systems, the setup time is related to the exchange (update) of the software transfer data/results, synchronization calculations, etc. Apart from the execution time of machine, the criteria, taking into account also other parameters (e.g. the number of used processors ) undoubtedly enable more accurate determination of the actual cost of operating of such a system.

Problem definition

The considered in the work problem of multi-machine tasks scheduling can be formulated as follows.

Problem: there is a set of tasks J = {1, 2, . . . , n}, given, which should be executed on identical machines (i.e. characterized with the same functional properties and equal capacities) from the set M = {1, 2, . . . , m}. Task i ∈ J requires at the same time execution of size i machines in

p i > 0 time units. Machines l, k ∈ M (l = k) are called neighbouring, if k = l +1 (l = 1, 2, ..., m-1) lub k = l -1 (l = 2, 3, ..., m)
. By s i,j we denote the time of machine setup after the completion of task i and before starting the task j (i, j ∈ J ), where both tasks are performed on the same machine. At the same time the following restrictions must be met: a) any machine, at any moment, cannot exercise more than one task, b) the task execution cannot be interrupted, c) each task is performed on the required number of neighboring machines, d ) between the tasks executed sequentially there must be setup of the machine performed.

We assume further that the number of machines is m max{size i : i ∈ J }.

The considered in the work problem (in short denoted by MPP) boils down to determining, for each task, a subset of machines and the starting moments of its execution (for each machine) satisfying the constraints (a)-(d), to optimize the adopted criterion. The first solution can be represented by a pair of Θ = (Q, S) is such that:

-Q = [Q 1 , Q 2 , ..., Q n ], where Q i (Q i ⊆ M, |Q i | = size i , ) is a set of machines
(assignment) on which the task i ∈ J , -S = [S 1,1 , ..., S 1,m , S 2,1 , ..., S 2,m , ..., S n,1 , ..., S n,m ], will be executed, where the element S i,j is the starting moment of task i execution on machine j (if j / ∈ Q i , then we assume S i,j = -∞).

By Ω let us denote the set of all solutions to MPP problem. For any solution

Θ ∈ Ω C max (Θ) = max{S i,j + p i : i ∈ J , j ∈ M}
is the moment of all task completion, and

M max = max{j : j ∈ Q i , i ∈ J }
the maximum number of machine from all allocated to tasks execution. Let

F (Θ) = C max (Θ) • M max (Θ). (1)
We consider the problem of scheduling of Textbf MPP tasks consisting in determination of solution Θ * ∈ Ω such that

F (Θ * ) = min{F (Θ) : Θ ∈ Ω}.
Example.There is a set of multi-machine tasks given J = {1, 2, 3, 4, 5, 6} to be executed on m = 4 machines from the set M = {1, 2, 3, 4}. The parameters of the individual tasks are shown in Table 1 (it is assumed that the setup times equal 0).

task (j) 1 2 3 4 5 6 pj 1 2 1 3 3 1 sizej 4 3 1 2 2 1 Figure 1 shows a Gantt diagram to certain solution Θ = (Q, S), where Q = [Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 ] and Q 1 = {1, 2, 3, 4}, Q 2 = {2, 3, 4}, Q 3 = {1}, Q 4 = {3, 4}, Q 5 = {1, 2}, Q 6 = {1}. The beginning times of execution of tasks S = (0, 0, 0, 0, -∞, 1, 1, 1, 1, -∞, -∞, -∞, -∞, -∞, 3, 3, 3, 3, -∞, -∞, 2, -∞, -∞, -∞). For this solution M max (Θ) = 4, C max (Θ) = 6,
therefore the value of the goal function is F (Θ) = 24. The work [START_REF] Lee | Scheduling one and two-processor tasks on two parallel processors[END_REF] proved that the problem of multi-machine task scheduling (without setup times) is NP-hard. Some lower and upper bound estimation of the criterion (1) were presented in work [START_REF] Drozdowski | Scheduling for parallel Processing[END_REF]. To solve the considered in the work MPP problem there was a simulated annealing algorithm presented. The starting moments of tasks execution S, at a fixed order, will be designated with the use of the algorithm solving two-dimensional packing problem (already onedimensional problem of packing is strongly NP-hard, Garey et all. [START_REF] Garey | Computers and intractability: a guide to the theory of NP-completeness[END_REF]).

Packing problem

Two-dimensional bin packing problem consists in such arrangement of rectangles, that the field occupied by them had minimum surface area. In the considered variant the rectangles, which contact the sides (hights) should be moved from one other to some extent. To put the matter more formally, it can be described as follows:

Problem: there is a set of rectangles R = {r 1 , r 2 , ..., r n } given. The sizes of the rectangle r k ∈ R are defined by a pair of numbers t k = (l k , w k ), where l k denotes the height, whereas w k -width of the rectangle. By d i,j we denote the distance by which the rectangle r i must be moved from r j , when they are placed right next to one another.

The considered problem is to find the G shape of a minimal surface area, within which it is possible to place all rectangles, i.e. minimization of the product,

Φ(G) = L • W, ( 2
)
where L is the height and W -the width of G. shape. In short, the considered in the paper packing problem will be denoted by BBP.

With each rectangle located within the G shape there is identified the position, in which it was enclosed, represented by the coordinates in a Cartesian coordinate system. The rectangles should be placed inside the G shape in such a way so as not to overlap one other, whereas their edges, representing respectively their length and width, were located parallely to the edges representing the length and width of G. shape. The interval in the horizontal plane, between any two rectangles r i , r j ∈ R cannot be less than the value d ri,rj . The example illustrating the above-described definitions and indications is given in Fig. 2. It is easy to observe the following relationship:

G X Y W L 𝑙 𝑘 𝑤 𝑘 𝑑 𝑖,𝑗
Corollary 1 The problem of multi-machine task scheduling MPP is equivalent to the problem of two-dimensional bin packing BBP.

Simulated annealing algorithm

The solution is represented in the form of a permutation (sequence) of rectangles

α = [r 1 , r 2 , . . . , r n ].
On its basis, using the packing strategy (the decoding solution), there is the value of the solution determined, i.e. the surface area of the G(α) shape in which rectangles are placed. Permutation Λ determines the order of rectangles in which they will be executed by the packing procedure. Let Ω be the set of all permutations of elements of R, i.e. the set of feasible solutions. The simulated annealing algorithm (SA) will be used to search through this set.

The key elements of the algorithm are: initial temperature, probability function, with which the solutions worse than the current one, are accepted and cooling scheme (decreasing of the value) of the temperature.

In each step of the algorithm there are some elements the solutions space Λneighborhood of the current solution considered. It is a set of elements, generated by small modifications of the considered solution. In the case where the current solution is worse than the one drawn from the neighborhood, it is replaced by it. With some probability, it is also possible to replace the current solution with the worse solution. This strategy is used not only to prevent stagnation but also to direct the trajectory of explorations in new areas of the solution space. The selection of temperature, cooling scheme and the probability function has influence on the frequency of accepting solutions worse than the current solutions, thus not only on the ability of the algorithm to leave local minima but also on the stability of the search itself.

Let α ∈ Ω be any starting solution. By N (α) we denote the neighborhood, whereas by t the control parameter (temperature).

Simulated annealing algorithm

α best ← α; repeat repeat Randomly determine element β from neighborhood N (α); if G(β) < G(α) then α ← β else if exp{-(G(β) -G(α))/t} > random then α ← β; if G(β) < G(α best ) then α best ← β; until Change control parameter ; Change control parameter t; until End Condition.
Neighborhood N (α) is generated by two types of random moves, namely -the swap and the twist type. Such neighborhoods are described in details in the works of Bożejko and Wodecki [START_REF] Bożejko | On the theoretical properties of swap multimoves[END_REF][START_REF] Bożejko | Solving Permutational Routing Problems by Population-Based Metaheuristics[END_REF][START_REF] Bożejko | Parallel genetic algorithm for the flow shop scheduling problem[END_REF]. The temperature, in turn, is changed according to the geometric distribution, i.e. t k+1 = t k • λ where lambda 0, 8 λ < 1.

Packing strategy

In the construction of rectangle packing procedures, i.e. in the process of determination of the shape (area) within which the rectangle is enclosed, there are different strategies applied. They differ from one another in the criteria for selecting the area (coordinates), in which the considered rectangle will be placed. In each of the used strategies we consider rectangles in the order they appear in the solution (permutation), at the same time trying to designate the space of the smallest possible surface area. A general description of the rectangle packing procedure is presented below.

Let α = [r 1 , r 2 , . . . , r n ] be a certain sequence (permutation) of rectangles. The rectangles are considered in the order they appear in α. In i-th iteration of the procedure (i = 1, 2, . . . , n) we will enclose the rectangle r i to G shape (area), in which there are already enclosed rectangles r 1 , r 2 , . . . , r i-1 . The set P i-1 = {z 1 , z 2 , . . . , z s }, where z i = (x i , y i ) contains positions on which the considered rectangle is inserted (we assume that at the starting moment the area G = ∅ and P 0 = {(0, 0)}). Next, rectangle r i is 'temporarily' placed into the G, area (shape) on each position from the set P i-1 . In this way, we generate a set of temporary areas O = {G 1 , G 2 , . . . , G s }. Depending on the used packing strategy, we choose from O corresponding area (i.e. the coordinates of the point where we insert a rectangler i ). If the area G l , was chosen, then: (i ) rectangle r i is placed into the G shape on the position z l = (x l , y l ), (ii ) the set of positions

P i = P i-1 \{z l }∪{(x l +d rc,ri •δ, y l ), (x l +d rc,ri •δ, y l +l i )},
where δ ∈ {0, 1} and takes the value 1, if the placed rectangle r i was moved by size d rc,ri , and position z l was created in the result of placing the rectangle r c . Otherwise δ = 0.

Next, we move to the next iteration. The procedure ends after placing all rectangles on their positions. When selecting position, in which the following rectangle is inserted, the following strategies were used:

1. Minimum area -the position, causing minimal enlargement of the area ( 2) is chosen, created as a result of the insertion of the considered rectangle, 2. Area-sizes -it is the development of the minimum area strategy. The value of the ratio is calculated, taking into account the surface area of the temporary (area) shape and the difference in the length of its sides. The purpose of application of this factor is to create areas having a uniform shape, 3. Roulette -the position is drawn in accordance with the principle of a roulette.

The probability of drawing is proportional to the size of the temporary area, 4. XYZ -similarly as in the case of a roulette but the drawing is carried out according to the uniform distribution.

In the simulated annealing algorithm, each of the above-described packing strategy was implemented. For each designated by the algorithm solution G there was packing accuracy calculated T /Φ(G), where T is the sum of all the areas of rectangles (the lower bound of function Φ).

Computational experiments

The computations were performed on a computer equipped with a 6-core Intel Core i7 X980 CPU (3.33GHz) and Linux operating system Ubuntu 12.04.5 LTS. For the purposes of computational experiments there were test data of varying degrees of heterogeneity randomly generated. The number of possible types of rectangles is 3(107), 5(101), 8(143), 10(134), 12(142), 15(129), 20(161). The total number of rectangles is given in brackets. The following values of the the SA algorithm parameters were adopted: initial temperature t 0 = 125, rate of temperature change λ = 0, 98, while the temperature change occurs every 100 iterations (return to the initial value every 500 iterations) and a maximum number of iterations is 5000. The results of the carried out study are presented in tables [START_REF] Błądek | On contiguous and noncontiguous parallel task scheduling[END_REF][START_REF] Bozoki | A branch-and-bound algorithm for thecontinous-process job-shop scheduling problem[END_REF][START_REF] Bożejko | On the theoretical properties of swap multimoves[END_REF][START_REF] Bożejko | Solving timetabling problems on GPU[END_REF]. Each row of the table contains the results obtained for one installation of test data: the number of types of rectangles (Types), the accuracy of the designated solution (Accuracy), the surface area of the shape (area) including the considered rectangles (Shape area), lengths of the shorter and the longer side of the surface area (Shorter, Longer side), the ratio of the longer side to the shorter side(Ratio) and the calculations time in seconds (Time). By far the best results for all instances of the problem were obtained by applying the strategy of minimum area (Table 2) in the packing procedure. For each example there was the shortest time of calculations observed. The parameter values Ratio were much greater than the value received after the application of other packing strategies (longer side was from 12.23 to 60.00 times larger than the shorter side). This strategy shows a clear preference of the surface area of the specified space. Using the Area-sizes ratio strategy (table 3) the results were only a little worse than those listed in the table 2. With the use of this strategy the smallest value of the Ratio coefficient for the number of types 3, 8, 10, 12, 15 and 20 were obtained. For the two types 3 and 5, the values were the same as in the strategy Roulette strategy. The value of the Ratio parameter close to 1.00 means that shape of the designated area is close to the square. The next two tables 4 and 5 include the results obtained after application of non-deterministic packing strategy in SW algorithm. The values of the solutions obtained with the use of the Roulette strategy (tab 4) were significantly worse than ones obtained with the use of Minimum area or area-size ratio strategy. Also, the time of compuatations was longer. For the two types 3 and 5 there were the biggest values of the Ratio parameter obtained. While applying the strategy XYP (the results -Table 5) there were the worst solutions, for the six types, obtained. This concerned both the computation time and the accuracy. Generally, deterministic strategies proved to be much better than probabilistic ones. If we take into account only the surface area, it is better to use the Minimum area strategy. On the other hand, when we are looking for an area of shape similar to a square, then it is better to use the Area-sizes strategy. For the determination of costs of multiprocessor tasks, the choice of packing strategy, depends on the relationship between the time of CPU utilization (computations time) and the number of processors used.

Simulated Annealing and Genetic Algorithm comparison

The canonical version of the Genetic algorithm was adopted. The algorithm has been configured as follows. For crossover, the PMX method was used, the Mutation was based on random chromosome fields exchange. The stopping condition was the maximum iteration number -5000. There are Accuracy, Ratio and Time results, for all strategies combined with Simulated Annealing and Genetic Algorithm compered in Table 6. The results are organized in four-row data structures, for all test data instances (3, 5, 8, 10 12, 15 and 20 types). Each row presents results for particular strategy and data set. The usage of Genetic algorithm in most cases improved the Accuracy results, especially for XYP strategy. Thirteen times the Ratio had greater value when Simulated Annealing algorithm was used, six times the result was equal. This means that structures builded with Genetic Algorithm had the shape close to square more often. The computation time was shorter in all cases for Area, Areasized and XYP strategy combined with Genetic Algorithm. Roulette packing strategy had shorter execution time with Simulated Annealing than with Genetic Algorithm.

Conclusions

In the paper there was the problem of multi-machine tasks scheduling with setup times and minimizing the product of the number of machines used and the length of time of all tasks execution considered. Designation of a solution can be brought down to a two-dimensional rectangles packing problem with minimization of the surface area of the occupied area (shape). The simulated annealing algorithm using different packing strategies was presented in the work, i.e. placing rectangles within the proposed area. The conducted computational experiments showed that the designated surface area (to solve the tasks scheduling problem) are only slightly greater than the sum of surface areas of all the rectangles (i.e. the lower bound).

Fig. 1 :

 1 Fig. 1: Gantt diagram for solution Θ.

Fig. 2 :

 2 Fig. 2: G shape with enclosed rectangles

Table 1 :

 1 Tasks parameters.

Table 2 :

 2 Minimum area strategy.

	Types Accuracy Shape area Shorter side Longer side	Ratio	Time
	3	0,920	489176	94	5204	55,36	56
	5	0.873	677100	150	4514	30.09	37
	8	0.850	766140	113	6780	60.00	117
	10	0.860	641200	229	2800	12.23	162
	12	0.842	689475	145	4755	32.79	171
	15	0.874	699470	226	3095	13.69	102
	20	0.841	892080	216	4130	19.12	160

Table 3 :

 3 Area-sizes ratio strategy.

	Types Accuracy Shape area Shorter side Longer side Ratio	Time
	3	0,861	522440	706	740	1,05	141
	5	0.851	694540	820	847	1.03	115
	8	0.832	782310	879	890	1.01	278
	10	0.803	686412	828	829	1.00	217
	12	0.816	711490	842	845	1.00	254
	15	0.819	746496	864	864	1.00	190
	20	0.804	933131	961	971	1.01	376

Table 4 :

 4 Roulette strategy.

	Types Accuracy Shape area Shortest side Longer side Ratio	Time
	3	0,220	2042350	1396	1463	1,05	142
	5	0.303	1947150	1379	1412	1.02	121
	8	0.183	3565740	1774	2010	1.13	317
	10	0.256	2151000	1434	1500	1.05	268
	12	0.193	3009260	1704	1766	1.04	315
	15	0.201	3043580	1689	1802	1.07	238
	20	0.148	5069220	2075	2443	1.18	444

Table 5 :

 5 XYP strategy.

	Types Accuracy Shape area Shorter side Longer side Ratio	Time
	3	0,134	3352340	1713	1957	1,14	180
	5	0.244	2423280	1366	1774	1.30	154
	8	0.157	4158810	1622	2564	1.58	423
	10	0.182	3035430	1305	2326	1.78	354
	12	0.167	3476110	1547	2247	1.45	417
	15	0.168	3632800	1353	2685	1.98	315
	20	0.151	4966610	1793	2770	1.54	604

Table 6 :

 6 Simulated Annealing and Genetic algorithms comparison

		Simulated Annealing	Genetic Algorithm
		Accuracy	Ratio	Time	Accuracy	Ratio	Time
	Area	0,920	55,36	55,7	0,920	55,36	15
	Area-sizes	0,861	1,05	140,9	0,894	1,05	16
	Roullette	0,220	1,05	141,7	0,256	1,16	149
	XYP	0,134	1,14	180,0	0,323	1,86	34
	Area	0,873	30,09	36,5	0,872	58,11	21
	Area-sizes	0,851	1,03	114,8	0,892	1,00	23
	Roullette	0,303	1,02	121,1	0,346	1,22	130
	XYP	0,244	1,30	154,3	0,462	1,24	37
	Area	0,850	60,00	117,4	0,898	15,31	30
	Area-sizes	0,832	1,01	278,1	0,866	1,00	63
	Roullette	0,183	1,13	317,4	0,236	1,44	345
	XYP	0,157	1,58	423,3	0,313	1,78	88
	Area	0,860	12,23	161,8	0,887	11,61	44
	Area-sizes	0,803	1,00	217,1	0,845	1,00	44
	Roullette	0,256	1,05	268,3	0,305	1,04	322
	XYP	0,182	1,78	354,2	0,418	1,22	69
	Area	0,842	32,79	170,9	0,903	20,30	54
	Area-sizes	0,816	1,00	254,4	0,852	1,00	52
	Roullette	0,193	1,04	315,2	0,261	1,04	340
	XYP	0,167	1,45	416,8	0,356	1,24	89
	Area	0,874	13,69	101,6	0,903	13,79	36
	Area-sizes	0,819	1,00	190,3	0,851	1,01	44
	Roullette	0,201	1,07	238,3	0,306	1,03	280
	XYP	0,168	1,98	315,0	0,390	1,33	60
	Area	0,841	19,12	160,1	0,853	9,58	54
	Area-sizes	0,804	1,01	376,1	0,835	1,01	91
	Roullette	0,148	1,18	443,5	0,218	1,07	489
	XYP	0,151	1,54	603,8	0,266	1,64	132


Micha L Piórek
email: michal.piorek@pwr.edu.pl

Mutual Information for Quaternion Time Series

Keywords: mutual information, quaternions, deterministic chaos, nonlinear time series analysis, human motion analysis, human gait data

Mutual Information method is a widely used method for estimation of time delay value in the process of time delay embedding. It's designed for a univariate scalar time series. In the real systems often many outputs of investigated system are available. In this case a multivariate time delay estimation method is necessary if one may require to perform the uniform time delay embedding. The special case of multivariate data is a kinematic time series(e.g quaternion time series or Euler angles time series). The main goal of this paper is to provide a method for this case: Mutual Information method's extension for quaternion time series. The results are also compared with previosly presented quaternion's angle method. The method was tested on the real kinematic data -the recordings of human gait.

Introduction

The time delay embedding is heavily explored area of nonlinear time series analysis and nonlinear dynamical systems areas. Usually the most common application of this method is when we want to discover the dynamics of underlying system from the univariate scalar time series created from the measurements of one of the investigated system's outputs. Taken's Embedding Theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF] implies, that one can reconstruct an equivalent dynamics from univariate time series using it's time delays. To carry out the embedding procedure, which should result in the reconstructed attractor in the output, two parameters need to be estimated: time delay and embedding dimension.

The time delay is a integer value describing which samples from the investigated time series we need to incorporate to time-lagged embedding vector -reconstructing the underlying phase space. There are few approaches of estimation of time delay value T d [START_REF] Chong-Zhao | Selection of embedding dimension and delay time in phase space reconstruction[END_REF]. One group of methods is series correlation approaches (autocorrelation, mutual information [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF] or high order correlations [START_REF]Using high-order correlations to define an embedding window[END_REF]). Second grop are approaches of phase space extension (fillfactor [START_REF] Buzug | Optimal delay time and embedding dimension for delay-time coordinates by analysis of the global static and local dynamical behaviour of strange attractors[END_REF], wavering product [START_REF] Buzug | Comparison of algorithms calculating optimal embedding parameters for delay time coordinates[END_REF] or average displacement [START_REF] Colins | Reconstruction expansion as a geometry-based framework for choosing proper delay times[END_REF]). There are available also multiple autocorrelation and non-bias multiple autocorrelation methods [START_REF] Huang | Selection of proper time-delay in phase space reconstruction of speech signals[END_REF]. The embedding dimension is an equivalent of the real underlying phase space dimension. It could be estimated using the false nearest neighbors method [START_REF] Brown | Determinig embedding dimension for phase-space reconstruction using a geometrical construction[END_REF] or it's extension -Cao's method [START_REF] Cao | Practical method for determining the minimum embedding dimension of a scalar time series[END_REF]. Another methods are also the saturation of system invariants method [START_REF] Abarbanel | Analysis of observed chaotic data[END_REF] or neural network approaches [START_REF] Maus | Neural network method for determining embedding dimension of a time series[END_REF].

The popularity of the univariate embedding may be caused by the fact, that according to the embedding theorem, for recovering dynamics only a univariate time series is needed. In fact, often many time series measured in the output of the test process are available. Since in multivariate case more data is available, it helps to establish more accurate embedding -in the sense of further predictions or in the presence of data noise. However, it brings a new dilemma: which quantities from multivariate time series to use and whether is better to use constant or non-constant embedding parameters for all quantities selected to embedding vector [START_REF] Cao | Dynamics from multivariate time series[END_REF].

The problem of multivariate time series embedding can be seen in terms of suitable conditioned embedding of the considered set of time series [START_REF] Vlachos | Nonuniform state-space reconstruction and coupling detection[END_REF]. In the related work there are two approaches of multivariate time series embedding: Uniform embedding and Non-uniform embedding. The uniform embedding scheme is more popular approach and assumes that embedding parameters: the time delay and embedding dimension are selected a priori and separately for each time series. The non-uniform embedding is based on the progressive selection of time delayed values from a set of candidate values(e.g X, Y , Z) and incorporation them to the embedding vector. In each step the most informative time delayed variables are chosen and then added to the time delay vector. As a selection criteria the mutual information between constructed embedding vector and the future state of the system is used [START_REF] Montalto | Mute: a matlab toolbox to compare established and novel estimators of the multivariate transfer entropy[END_REF] [START_REF] Vlachos | Nonuniform state-space reconstruction and coupling detection[END_REF].

A particular case of multivariate time series is a rotational data time series. There are three main parametrization of rotations: matrix of rotation, Euler angles and Quaternions. Basing on this fact one may record and construct rotational time series according to the one of the above parameterizations. In this paper quaternion rotational time series is considered.

The main goal of this work is to propose time delay estimation method for uniform time delay embedding of multivariate quaternion rotational data. The proposed approach bases on mutual information approach and it's re-designed for quaternion kinematic time series. The presented method could be used in the further time delay embedding and nonlinear analysis aimed to detect deterministic chaos properties in the investigated data. The author would like to underline that the considered method allows to estimate the time delay value staying in quaternion domain, which should help to keep physical sense of the kinematic data.

The paper is organized as following: in the second section the applicability of mutual information method for quaternion data is discussed and proposed approach is described. It also includes the information about investigated quaternion time series and how K-Means algorithm is applied. The third section presents the numerical results. The conclusions are presented in section four.

Mutual Information approach for Quaternions

Quaternions are computationally efficient parametrization of rotational data. They are an extension of complex numbers defined as following:

q = [w, (x, y, z)] = w + ix + jy + kz (1)
where: w represents a real part and v = (x, y, z) is called a vector part (i, j and k are equivalents of imaginary unit).

The details of Quaternions algebra widely used in the parametrization of rotations is well described in the related work(e.g [START_REF] Jablonski | Quaternion dynamic time warping[END_REF]). In the scope of our interests are unit quaternions which describe the rotation in 3D space:

∥q∥ = 1 (2)
Where quaternion norm is defined by:

∥q∥ = √ w 2 + x 2 + y 2 + z 2 (3)
We assume that the method is designed for the following quaternion time series formed by unit quaternions as following:

Q(n) = (q 1 , q 2 , ..., q N ) = (w 1 + ix 1 + jy 1 + kz 1 , ..., w N + ix N + jy N + kz N ) (4)

Mutual Information -existing approach

The mutual information is a measure which describes the general dependence of two variables. The definition comes from Shannon's information theory, which gives the formalism of measuring information spreading. Frasser proposed to use this approach in time delay estimation process [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF].

Let's assume that there are two nonlinear systems: A and B. The outputs of these systems are denoted as a and b, while the values of these outputs are represented by a i and b k . The mutual information factor describes how many bits of b k could be predicted where a i is known.

I AB (a i , b k ) = log 2 ( P AB (a i , b k ) P A (a i )P B (b k ) ) , (5)
where P A (a i ) is the probability that a = a i and

P B (b k ) is the probability that b = b k and P AB (a i , b k ) is the joint probability that a = a i and b = b k .
The average mutual information factor can be described by:

I AB (T ) = ∑ ai,b k P AB (a i , b k )I AB (a i , b k ). (6)
In order to use this method to assess the correlation between different samples in the same time series, the Average mutual information factor is finally described by the equation:

I(T ) = ∑ N n=1 P (S(n), S(n + T )) log 2 ( P (S(n),S(n+T )) P (S(n))P (S(n+T )) ) . ( 7
)
Fraser and Swinney ( [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF]) propose that T m where the first minimum of I(T ) occurs as a useful selection of time lag T d . This selection guarantees that the measurements are somewhat independent, but not statistically independent. In case of absence of the average mutual information clear minimum, this criterion needs to be replaced by choosing T d as the time for which the average mutual information reaches four-fifths of its initial value:

I(T d ) I(0) ≈ 4 5 . ( 8
)

Mutual Information extension for quaternion time series

The average mutual information method for univariate time series consists of 2-dimensional adaptive histogram and that is the problem in it's application to multivariate(quaternion's case). The empirical histogram is straightforward to estimate for univariate time series, however for quaternion's time series it's not trivial. Computation of multivariate histogram is exhaustive process and in the result one may obtain the histogram empty in some places.

In the current approach instead of multidimensional histogram for quaternions we propose here to use histogram based on clusters. The whole quaternion time series is initially clustered into k-groups (where k is defined a priori). The obtained clusters are treated as an equivalent of histogram bins. Further in empirical histogram estimation, instead of computing the probability of belonging to the histogram's bins, the probability of belonging to the clusters is being computed. C ← Initialize clusters randomly 5:

Algorithm 1 k Means algorithm

i ← 0Number of iterations 6:

while Ci! = Ci+1 do 7:

Assign instances to the closest cluster center 8:

Update cluster centers based on the assignment 9:

end while 10: end procedure Data clustering as a part of machine learning and data sciences is an actively investigated field of science. There are many available clustering techniques. The review of clustering methods is presented in the related work e.g [START_REF] Rokach | Clustering methods[END_REF]. In the presented approach K-means clustering algorithm was selected as a the simplest and commonly used method. The main goal of this work is not to examine the efficiency of clustering approaches but to provide the mutual information estimation technique used clusters based histogram. In the oder hand, the author see the underlying potential in investigation of the impact of clustering method selection on the general algorithm's performance.

The K-means algorithm was described by MacQueen [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. It partitions the data into K clusters (C 1 , C 2 , ..., C k ), represented by their centers. The center of each cluster, until converge, is calculated as the mean of all the samples belonging to that cluster. Initially the centers are selected randomly. Then, in each iteration each sample is assigned to the closest cluster center according to the Euclidean distance. Then the centers are re-calculated. The whole procedure is repeated until the convergence criteria is fulfilled (e.g there is no relocation of the centers in new iteration) [START_REF] Rokach | Clustering methods[END_REF]. The whole procedure is presented in the pseudo-code 1. Finally, the mutual information algorithm for a quaternion time series is described by the pseudo code [START_REF]Using high-order correlations to define an embedding window[END_REF]. Initially it treats quaternion time series as a 4-dimensional time series and partitions it into K-clusters. The next step is an estimation of probability of belonging samples to the each cluster which is an equivalent of estimation of probability of belonging samples to histogram bins in a standard version of the algorithm.

Algorithm 2 Mutual information algorithm for quaternion time series

The author sees the potential advance in the method using an algorithm of clustering which partitions quaternion time series staying in quaternion domain. It will be the subject of further research. It is also worth to investigate how the number of the clusters impacts the performance of the algorithm.

Numerical results

The method was tested on live kinematic data recorded in the Human Motion Laboratory (HML) of the Polish-Japanese Institute of Information Technology. The recordings were performed using the Vicon Motion Kinematics Acquisition and Analysis System equipped with 10 Near InfraRed cameras. The cameras were attached to a suit which was worn by a subject. Gait sequences were recorded in Euler angles and then converted to quaternions. Six kinds of time series were recorded -movements of femurs, tibias and feet (left and right). The method was tested on the data recorded from the treadmill walking of two healthy patients. The number of clusters in K-means algorithm was set to 7. The designed method was additionally compared with the quaternion angle method presented by the author in the same conference last year [START_REF] Piorek | Computer Information Systems and Industrial Management: 14th IFIP TC 8 International Conference, CISIM 2015[END_REF]. The results for the designed method are presented using solid line, where the results performed using quaternion angle method are presented using dashed line. Constantly the first local minima of the mutual information functions(time delay selection criteria) was marked by vertical lines. All estimated time delays are gathered in the table1.

Conclusion

The main goal of this article was to present a time delay estimation method for a quaternion time series. The approach extends the existing mutual information approach for quaternion time series by incorporation of K-means clustering for multivariate data instead of the empirical histogram. The method might be a first step to perform time delay embedding staying only in the quaternion domain, which will be a field of author's further research. From visual inspection one can see that the results from the proposed method are in the same range as the results coming from previously investigated quaternion angle method. It is worth to underline that the differences in the result are expected, since quaternion angle method bases only on a part of quaternion's information where the new method utilize the whole information carried by a quaternion.

The field of further interests should be also the impact of clustering parameters on the methods performance and the analysis of the quality of embedding using the proposed approach.

Fig. 1 .

 1 Fig. 1. Mutual information dependency for the patient A

Fig. 2 .

 2 Fig. 2. Mutual information dependency for the patient B

Table 1 .

 1 Time delay estimation comparison for subject A and subject B

	body part	Td MI-angle(Patient A)	Td MI-quat(Patient A)	Td MI-angle(Patient B)	Td MI-quat(Patient B)
	left femur	25	28	13	6
	right femur	24	18	11	5
	left foot	17	17	14	8
	right foot	30	32	13	8
	left tibia	38	26	24	22
	right tibia	25	25	24	22

Acknowledgments. This work has been supported by the "M loda Kadra" project from the Wroc law University of Science and Technology.

The author also would like to acknowledge the Polish-Japanese Institute of Information Technology for the gait sequences recorded in the Human Motion Laboratory (HML).