extfieldsofstudy:
Mathematics
Computer Science

text:
An Efficient Meshfreee Implicit Filter for Nonlinear Filtering Problems

In this paper, we propose a meshfree approximation method for the implicit filter developed in [2], which is a novel numerical algorithm for nonlinear filtering problems. The implicit filter approximates conditional distributions in the optimal filter over a deterministic state space grid and is developed from samples of the current state obtained by solving the state equation implicitly. The purpose of the meshfree approximation is to improve the efficiency of the implicit filter in moderately high-dimensional problems. The construction of the algorithm includes generation of random state space points and a meshfree interpolation method. Numerical experiments show the effectiveness and efficiency of our algorithm.

where Φ(·) is a test function, and Z k denotes the space of all Y 1:k -measurable and square integrable random variables.
Tremendous efforts have been made to solve nonlinear filtering problems in the last few decades. Two of the well-known Bayesian filters are extended Kalman filters (EKFs) [3,10,16,17,19], and particle filters [5,7,8,13]. The key ingredient of the EKFs is the linearization of both f and g in (1.1), so that the standard Kalman filter can be applied directly. Thus, if the nonlinearity of the state and the observation systems is not severe, then the EKFs can provide efficient and reasonable inferences about the state, otherwise, the performance of the EKFs can be very poor. For particle filters, the central theme is to approximate the desired posterior PDF of the state by the empirical distribution of a set of adaptively selected random samples (referred to as "particles"). The particle filter method is essentially a sequential Monte Carlo approach, which requires no assumption on the linearity of the underlying system. As such, with sufficiently large number of particles, it is capable of providing an accurate approximation of the posterior PDF for a highly nonlinear filtering problems. However, there are some fundamental issues concerning the efficiency and robustness of particle filters [9]. For example, since the empirical PDF is constructed based on particles with equal weights after resampling, the particle filter still needs a lot of samples in order to accurately approximate the target distribution.
To overcome such a disadvantage, the authors proposed a new nonlinear filter named "implicit filter" [2] . This approach adopts the framework of Bayesian filtering, which has two stages at each time step, i.e., prediction and update. At the prediction stage, we estimate the prior PDF p(X k |Y 1:k−1 ) of the future state X k given the current available observation information Y 1:k−1 ; at the update stage, we update the prior PDF by assimilating the newly received data Y k to obtain the estimate of the posterior PDF p(X k |Y 1:k ). The implicit filter is distinguished from the particle filters by the use of interpolatory approximations to the prior and posterior PDFs. Specifically, in the particle filter, p(X k |Y 1:k−1 ) is approximated by explicitly propagating the samples of the current state X k−1 |Y 1:k−1 through the nonlinear state equation X k = f k (X k−1 , w k−1 ), and constructing the empirical PDF of X k |Y 1:k−1 . In the implicit filter, the interpolation of p(X k |Y 1:k−1 ) requires its function values at a set of grid points of the future state X k . Under the condition that X k = x ∈ R d , we solve implicitly the state equation x = f k (X k−1 , w k−1 ) given a set of Monte Carlo samples of w k−1 , so that the value of p(X k = x|Y 1:k−1 ), at the grid point of x, can be estimated by averaging the function values of p(X k−1 |Y 1:k−1 ) at all the solutions of the state equation. As an implicit scheme, the implicit filter has a stabilizing effect which provides more accurate numerical approximations to the solution of the nonlinear filtering problem than the particle filter method [2].
The main challenge of the implicit filter method is that the conditional PDF of the nonlinear filtering solution is estimated at grid points. As such the method suffers the so called "the curse of dimensionality" when the dimension of the state variable is high. In addition, the efficiency of the method may be significantly reduced when the domain of the PDF is unbounded. In this paper, we propose to construct a meshfree implicit filter algorithm to alleviate the aforementioned challenges. Motivated by the particle filter method, we first generate a set of random particles and propagate these particles through the system model and use these particles to replace the grid points in the state space. After that we generate other necessary points through the Shepard's method which constructs the interpolant by the weighted average of the values on state points [12]. In order to prevent particle degeneracy in the generation of random state points, we introduce a resample step in the particle propagation. In addition we choose state points according to the system state, which make them adaptively located in the high probability region of the PDF of state. In this way, we solve the nonlinear filtering problem in a relatively small region in the state space at each time step and approximate the solution on a set of meshfree state points distributed adaptively to the desired PDF of the state. Furthermore, since we approximate the PDF as a function on each state point, instead of using state points themselves to describe the empirical distribution, the implicit filter algorithm requires much fewer points than the particle filter method to depict the PDF of the state.
The rest of this paper is organized as follows. In §2, we introduce the mathematical framework of the Bayesian optimal filter. In §3, we construct meshfree implicit algorithm. In §4, we demonstrate the efficiency and accuracy of our algorithm through numerical experiments. Finally, §5 contains conclusions and directions for the future research.
2. Bayesian optimal filter. For m, n ∈ N + , let X m:n and Y m:n denote the σ fields generated by {X m , X m+1 , . . . , X n } and {Y m , Y m+1 , . . . , Y n }, respectively. For k = N + , we use x k to represent a realization of the random variable X k , and define for notational simplicity. It is easy to see that the dynamical model in (1.2) is Markovian in the sense that We also know that the measurements are conditionally independent given x k , i.e., The Bayesian optimal filter constructs the conditional distribution p(x k |Y 1:k ) recursively in two stages: prediction stage and update stage.
For k = 1, 2, · · · , assume that p(x k−1 |Y 1:k−1 ) is given. In the prediction stage p(x k |Y 1:k−1 ) is evaluated through the Chapman-Kolmogorov formula: (2.1) In the update stage, the prior PDF obtained in (2.1) is used to obtain the posterior PDF p(x k |Y 1:k ) via the Bayes' formula: 3. The meshfree implicit filter. In this section, we construct the meshfree implicit filter algorithm. The algorithm is based the implicit filter algorithm on grid points [2]. The implicit filter algorithm introduced in [2] is developed from the general framework of the Bayesian optimal filter discussed above, in which the primary computational challenge is the numerical approximation of the term p(x k |x k−1 ) in (2.1).
3.1. The prediction stage. For k = 1, 2, · · · , the goal of this stage is to approximate the prior distribution p(x k |Y 1:k−1 ) of the state X k , given the posterior distribution p(x k−1 |Y 1:k−1 ) of the state X k−1 . Due to the the fact that the prior PDF p(x k |Y 1:k−1 ) derived in identity (2.1) can be rewritten as represents the expectation with respect to the white noise w k−1 , and the PDF p(x k |x k−1 , w k−1 ) is 2) can be viewed as a generalization of the Dirac delta function in the space R d × R d × R r , where the mass is located according to the state equation Note that the estimation of (3.1) requires an approximation to the expectation of the white noise w k−1 , and define an approximation to p(x k |x k−1 , w k−1 ) as To construct an interpolation of p(x k |Y 1:k−1 ), the next step is to approximate is the approximation error. Then, by further fixing w k−1 = w j k−1 , the location of the mass can be obtained by implicitly solving the state equation which is the reason we named the approach the implicit filter. Now substituting x i k into (3.4), and using the same sample set {w j k−1 } M j=1 as above to approximate the integral on the right hand side of (3.4), we obtain , we obtain the following iterative numerical scheme for constructing an approximation, denoted by (x i k |Y 1:k−1 ), of the prior PDF p(x i k |Y 1:k−1 ), i.e., In our previous work [2], the subsets H k , for k = 0, 1, . . ., were defined by a full tensor product mesh, denoted by It is simple to implement but has several significant disadvantages. First, at each time step, one needs to approximate the prior PDF p(x k |Y 1:k−1 ) at a total of N (1) × · · · × N (d) grid points which grows exponentially as the dimension d increases. This is also known as "the curse of dimensionality". On the other hand, since the construction of M is not informed by the target PDF, the domain [a 1 , b 1 ] × · · · × [a d , b d ] needs to be defined sufficiently large, so as to capture the statistically significant region of the PDF. This may lead to a great waste of computation effort in the low probability region of p(x k |Y 1:k−1 ).
To alleviate such disadvantages, we propose to develop a distribution-informed meshfree interpolation approach to efficiently approximate the prior PDF. The central idea of the generation of random points for the state variable is to build a set of points, denoted by H k , according to the state distribution. To begin with, we generate of N random samples from the initial PDF p 0 of the initial state: If the initial PDF p 0 is close to the true state distribution, it's obvious that our random state points are more concentrated near the target state. For k = 1, 2, · · · , K, we propagate points and approximate the conditional PDF p(x k |Y 1:k−1 ) on H k with the scheme given by (3.8). In this way, the random points in H k move according to the state model. As opposed to particle filter methods, which use the number of particles to represent empirical distributions and require a large number of particles to follow the state distribution, in the implicit filter method we provide an approximation of the value of the PDF at each state point. Therefore, much fewer points are needed to describe the state PDF and the random state points are not necessary to accurately follow the state distribution.
3.2. The update stage. By incorporating the new data Y k , we update the prior PDF p(x k |Y 1:k−1 ) at each grid point x i k , using the Bayesian formula, to obtain where (x i k |Y 1:k−1 ) is given in (3.8), C k is the normalization factor, and R i k|k : is the approximation error. By neglecting the error term R i k|k in (3.10), we obtain the following iterative numerical scheme for the update stage on H k , i.e., where (x i k |Y 1:k ) is desired the approximation of the posterior PDF p(x i k |Y 1:k ). Next, we use interpolation methods to construct the approximation (x k |Y 1:k ) of is the set of basis functions. Since the state points in H k are generated randomly in the meshfree framework, standard polynomial interpolation [2] is unstable due to the uncontrollable Lebesgue constant. Instead, we propose to use the Shepard's method [12], which is an efficient meshfree interpolation technique, to construct the interpolant (x k |Y 1:k ). The basic idea of the Shepard's method is to use the weighted average of { (x i k |Y 1:k )} x i k ∈H k in the interpolating approximation. Specifically, for a given point x k ∈ R d , we re-order the points in H k by the distances to x k to get a sequence {x where · is the Euclidean norm in R d . Then, for a pre-chosen integer L we use the where the weight h l (x k ) is defined by is the error of the Shepard's interpolation. We assume that p(x k |Y 1:k ) has bounded first order derivative. For each pair p(x k |Y 1:k ) and p(x k . It is reasonable to assume that in high probability region of the derivative p (z|Y 1:k ) is large. It's worth pointing out that the random state points generated in this algorithm are concentrated in the high probability region. Thus, if x k lies in the high probability region, the distance is small, which balances the error brought by the large derivative. On the other hand, if x k lies in the low probability region, although the distance x k − x (l) k is relatively large, the approximation error (3.14) is still small due to the small value of the derivative p (z|Y 1:k ).

3.3.
Resampling. Similar to the particle filter method, the above random state points generation suffers from the degeneracy problem for long term simulations, especially for high-dimensional problems. After several time steps, the probability density tends to concentrate on a few points which dramatically reduces the number of effective sample points in H k .
In this work, we propose an occasional resampling procedure to address these problems and rejuvenate the random points cloud. At the time step k − 1, the resampling procedure takes place after we obtain (x k−1 |Y 1:k−1 ), in order to remove the degenerated points in H k−1 using the information provided by (x k−1 |Y 1:k−1 ). Specifically, the first step is to develop a degeneracy metric to determine the necessity of doing resampling. To this end, we define the following degenerated subset S k−1 ⊂ H k−1 , where ε > 0 is a user-defined threshold. We also define J (S k−1 ) := {i = 1, . . . , N |x i k−1 ∈ S k−1 } to be the index set of S k−1 . Then, the degeneracy of H k−1 can be measured by the ratio #(S k−1 )/#(H k−1 ) ∈ [0, 1], where #(·) denotes the number of points in a set. If the ratio is smaller than a threshold τ ∈ [0, 1], then we will skip the resampling step and propagate H k−1 to get H k ; otherwise, the set H k−1 is considered degenerated, and the resampling procedure is needed.
In resampling, instead of propagating H k−1 to H k , we aim at constructing an and propagate H k− 1 2 through the state model (1.2) to obtain H k . According to the definition of S k−1 in (3.15), we consider the state points in H k−1 \S k−1 are in the statistically significant region of (x k−1 |Y 1:k−1 ), so that we first put those points in H k− 1 2 , i.e., For the state points in S k−1 , we replace them by generating new samples from (x k−1 |Y 1:k−1 ) using the importance sampling [6], i.e., As a result, the resampling procedure helps us remove the state points with low probabilities, and makes the state point set H k concentrated in the high probability region of the posterior PDF (x k−1 |Y 1:k−1 ) at each time step.
3.4. Summary of the algorithm. Finally, we summarize the entire meshfree implicit filter algorithm introduced in §3.1- §3.3 in Algorithm 1 below.

Numerical experiments.
In this section, we present two numerical examples to examine the performance of our meshfree implicit filter method. In Example 1, we use a two dimensional nonlinear filtering problem to show the distributions of the random points H k . In Example 2, we solve a three dimensional bearing-only tracking problem, which is a six dimensional nonlinear filtering problem. For this higher dimensional problem, we compare the accuracy and efficiency of our meshfree implicit filter method with the extended Kalman filter and the particle filter. Example 1. In this example, we consider the two dimensional noise perturbed tumoral growth model [21] where W t is a two dimensional standard Brownian motion and σ = (0.01, 0.01) T . The state process X t = (X 1 t , X 2 t ) T is a two dimensional vector, F (X t ) := (f 1 (X t ), f 2 (X t )) T is defined as .
Here, f 1 models the Gompertzian growth rate of the tumor and f 2 gives the degree of vascularization of the tumor which is also called " angiogenic capacity". To approximate the state variables, we discretize the dynamic system (4.1) in time and obtain a discrete state model Here, ω k is a two dimensional zero mean Gaussian white noise process with covariance Q = I∆, where I is the 2 × 2 identity matrix and ∆ is the time partition stepsize. The measurement of the state model is given by where v k is a two dimensional zero mean Gaussian white noise process with covariance Λ = I∆, I is a 2 × 2 identity matrix and R = (0.1, 0.1) T . In the numerical experiment, we use uniform time partition with stepsize ∆ = 0.2 and simulate the state process for K = 40 with initial state X 0 = (0.8, 0.3) T and parameters α 1 = 1, α 2 = 0.2, α 3 = 0.2. At time step k = 0, we initialize the prior PDF p 0 by N (X 0 , Σ), whereX 0 = (0.78, 0.32) T and In Figure 4.1, we plot 1500 random samples generated from the initial PDF p 0 , which are our initial random points H 0 . To present the accuracy of the algorithm, we show the simulation of the tumoral growth states in Figure 4.3. The black curves are the true X 1 and X 2 coordinate values of the tumoral growth states, respectively. The blue curves show the simulated states obtained by using the meshfree implicit filter method.
Example 2. In this example, we study a six dimensional target tracking problem. In Figure 4.4, the target, denoted by the red line, moves in the three dimensional space and two platforms on the ground, denoted by pentagons, take angular observations of the moving target.
The state process X k = (X 1 k , X 2 k , X 3 k , X 4 k , X 5 k , X 6 k ) T is described by the following dynamic model where (X 1 , X 2 , X 3 ) describes the position of the moving target which is controlled by parameters (X 4 , X 5 , X 6 ). The system noise ω k = (ω 1 k , ω 2 k , ω 3 k , ω 4 k , ω 5 k , ω 6 k ) T is a zero mean Gaussian white noise process with covariance Q . = I∆, I is the 6 × 6 identity matrix and ∆ is a given time period, σ = (0.1, 0.1, 0.1, 0.01, 0.01, 0.01) T is a constant vector and f is given by where v k is a 4 dimensional zero mean Gaussian white noise process with covariance Λ = I∆, I is a 4 × 4 identity matrix, R = (0.6, 0.6, 0.6, 0.6) T , (a 1 , b 1 ) and (a 2 , b 2 ) are locations of two observers. We choose ∆ = 0.3, α = 3, v 1 = v 2 = v 3 = 0.05 . Also, we assume that platforms are located at (a 1 , b 1 ) = (16, 6), (a 2 , b 2 ) = (8, 15) and the initial sate is given by a Gaussian N (X 0 , Σ) where X 0 = (2, 2, 1, 0.4, 0.4, 0) T and The target will be observed over the time period 0 ≤ t ≤ 15. In the numerical experiments, we compare the performance of our meshfree implicit filter with the extended Kalman filter and the particle filter. In particular, we compare the estimated mean values of the states process along each dimension in Figure 4.5. In the particle filter method, we choose 15, 000 particles. In the meshfree implicit filter method, we choose the number of state points to be N = 4, 000 and the number of random samples in the implicit filter Monte Carlo simulation to be M = 6. The black curves in Figure 4.5 show the real states process along each direction, the green curves give the estimated means obtained by the extended Kalman filter method, the red curves give the estimated means obtained by the particle filter method, and the blue curves give the estimated means obtained by the meshfree implicit filter. We also plot the L 2 error err k corresponding to all three methods in figure 4.6. As we can see from figure 4.5 and 4.6, the implicit filter and the particle filter are much more accurate than the extended Kalman filter and the implicit filter is the most accurate approximation in this experiment.
To further compare the efficiency between the meshfree implicit filter and the particle filter, we repeat the above experiment over 50 realizations and show the average CPU time and the corresponding global root mean square error err G defined by where err k (j) is the L 2 error of the j-realization at time step k. In table 4.1, we can see that with 15, 000 particles, the CPU time of the particle filter method is comparable to that of the implicit filter with 4, 000 random state points, but the global RMSE of the particle filter is more than doubled the RMSE of the implicit filter. With 25, 000 particles, the particle filter method achieves an accuracy comparable to the implicit filter, but at a significantly higher cost.

Conclusions.
In this work, we proposed an efficient meshfree implicit filter algorithm by evaluating the conditional PDF on meshfree points in the state space. These meshfree points are chosen adaptively according to the system state evolution. We also apply Shepard's method as the meshfree interpolation method to compute interplants with random state points. In order to address the degeneracy of the random points, we use importance sampling method to construct a resample step. Numerical examples demonstrate the effectiveness and efficiency of our algorithm. In   the future, we plan to perform a rigorous numerical analysis for the meshfree implicit filter algorithm.


extfieldsofstudy:
Computer Science
Medicine

text:
Multiple Human Tracking Using Binary Infrared Sensors

To create a context-aware environment, human locations and movement paths must be considered. In this paper, we propose an algorithm that tracks human movement paths using only binary sensed data obtained by infrared (IR) sensors attached to the ceiling of a room. Our algorithm can estimate multiple human movement paths without a priori knowledge of the number of humans in the room. By repeating predictions and estimations of human positions and links from the previous human positions to the estimated ones at each time period, human movement paths can be estimated. Simulation-based evaluation results show that our algorithm can dynamically trace human movement paths.

Introduction
Human tracking technologies have attracted considerable attention over the years. For example, with real-time human tracking technology, air conditioners and lights can be controlled smartly by considering the requirements of each human. In addition, elderly people or children can be monitored for safety reasons [1,2]. Consequently, many human tracking systems and algorithms have been proposed [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. One popular method involves camera-based systems [3][4][5][6][7][8][9][10]. However, they are often not acceptable for monitoring elderly people or other similar applications whose objective is to observe humans without invading their privacy. This is because people may be uncomfortable with monitoring using video cameras, even if the original image data are modified and not used directly to perform human tracking. In addition, the size of video data is large compared to other sensor data; therefore, OPEN ACCESS camera-based systems are not suitable for long-term observations. In other popular methods, IC tags or radio frequency identification devices (RFIDs) [11,12] are used. A person wears a RFID, which is then detected by nearby RFID receivers located on the floor or ceiling of the monitored area. Thus, using data stored in the RFID receivers, the person can be tracked. Although the RFID-based systems can be considered to enhance privacy, they force people to carry the RFID devices.
For these reasons, some human tracking methods use infrared (IR) sensors [13][14][15][16][17][18][19][20]. The IR sensor produces a "1" if it detects a human, and a "0" otherwise [21]. Thus, using the IR sensors, the movement paths of subjects can be estimated without sacrificing their privacy. In addition, IR sensors are low-cost and their installation is relatively simple. However, IR sensors have a deficiency-they react regardless of whether they detect one or multiple persons. Therefore, almost all human tracking methods using IR sensors assume that the number of humans in the room is known beforehand. However, from the perspective of practical application development, the number of humans in a room is often unknown.
We previously proposed an algorithm that can simultaneously estimate both the number of humans and their movement paths [22,23], however, the method still had some drawbacks. Thus, the algorithm is invoked after all sensed data are collected, making it difficult to estimate the human movement in real time. In addition, the estimated human position must be the position of a fired IR sensor or the midpoint of the overlapped detection areas of multiple fired IR sensors; thus, the method lacks flexibility. To address these issues we now propose a novel method that can overcome the aforementioned problems [24]. The method dynamically estimates human positions using the weighted centers of grouped fired IR sensors, instead of simply using the position of a fired IR sensor or the midpoint of the overlapped detection areas of multiple fired IR sensors in the previous methods. Thus, by using the new method, we can estimate multiple human movement paths in a timely manner, and the estimated human positions are not restricted by the locations of the IR sensors, thus improving the human tracking accuracy.
The remainder of this paper is organized as follows: in Section 2, a model of the assumed IR sensor system is described. Next, the details of the proposed algorithm are explained in Section 3. Then, evaluation results are shown in Section 4. Finally, our conclusions and future studies are described in Section 5.

Model of IR Sensor System
In this paper, we consider a system constructed using IR sensors and a PC. The IR sensors are randomly installed on the ceiling of a room, but their locations are known. In addition, each IR sensor has a non-directional and circular detection range, the radius of which is r. Here, r is typically several meters, if we use commercially available IR sensors [21]. We assume r = 2.0 m in this paper. The sensor outputs binary data, i.e., it produces a "1" if it detects one or more humans and a "0" otherwise. Moreover, all the IR sensors are connected to the PC using a wired or wireless communication network. Thus, data obtained from the IR sensors are collected in the PC. The data sampling rate is 6 Hz, and the sensed data are saved on the PC with a time stamp. In addition, the positional coordinates of each sensor is known. Furthermore, the detection ranges of the IR sensors should cover the entire area to detect human movements. Note that some IR sensor arrays have been proposed recently [25], and they can detect human movements, but their detection range is typically several meters. Thus, even if using the IR sensor array, it is difficult to cover a large area. Furthermore, the method proposed in this paper can naturally handle the IR sensor array by treating it as a set of densely deployed IR sensors, if the sensor array produces a set of {0, 1}-data depending on human detection.

Proposed Algorithm
We can trace a human movement path by connecting the positions of the fired sensors step by step, if the position coordinates of each sensor is known. However, this method presents the problem that if some of the IR sensors are close together, they often detect the same event and output "1" data at the same time, which could make it difficult to track the human movements using the simple method mentioned above. This is especially true when there is more than one person in the room. To solve this problem, we propose a heuristic algorithm that can estimate multiple human movement paths using only the binary sensed data.

Notations
First, to explain the procedures of our algorithm, we define the variables and parameters as follows: · |*|: The number of elements in list *. · CC(t): A set of the weighted center coordinates of the clusters at time t. CC(t) = {cci(t)|i = 1, …, Cmax(t)}, where cci(t) = (x, y) and they are referred to as cci(t).x and cci(t).y, respectively. · td.Path: A list of the estimated route coordinates of target human "td.ID." By connecting all elements in this list from the first to the last, the estimated movement path of human "td.ID" is obtained. · td.TTL: A lifetime of target human "td.ID". This is an integer value. The initial value of td.TTL is "1," and its maximum value is TTLMAX. · TDL: A list of target humans who currently exist in the room. · TTLMAX: A constant integer value. The maximum value of TTL (time to live). TTLMAX > 0. · WL(t): A weight set of the sensors at time t. WL(t) ={ws(t)|s = 1, …, S}. · ws(t): A weight of sensor s at time t. ws(t) = ∑ d i |s = 1, … , S. This is used to calculate the weighted center coordinate of a cluster, to which sensor s belongs. The coordinate calculation method will be defined later. · WS1: Window size for the number of sensed raw datasets to which the logical-OR operation is applied. This is a given integer value. WS1 > 0. · WS2: Window size or the number of elements (coordinates) that construct movement vector mvt(td.ID). This is a given integer value. WS2 > 1.

// Location estimation step
Apply logical-OR operation and summation to RDL, and obtain ORD(t) and WL(t); call Clustering(ORD(t)) shown in Figure 3 The main procedure of the proposed algorithm is shown in Figure 1. It mainly consists of two processing steps: the location estimation step and path estimation step. The former is described as a part of the main procedure, whereas the latter step will be shown later. At each time t, the two steps are invoked sequentially, and a link from the current to the next human position is estimated. First, as the location estimation step, the candidates for human positions are listed using a clustering technique. Then, in the path estimation step, the links are determined by selecting reasonable positions among the candidates obtained in the location estimation step. At time t, the position of each human is estimated using the candidate sets, i.e., the sets of the weighted centers of clusters, obtained from t-WS2 + 1 to t. Here, WS2 is relatively small. Thus, our algorithm can trace the human movement paths almost in real time, without knowing the number of humans in the room. Our previous algorithm could not do this [22,23]. In addition, the framework of the proposed algorithm is the same as that of the algorithm introduced in [24], except that both the location estimation and path estimation steps are modified to improve human tracking accuracy. In the following subsections, the detailed procedure of each step will be explained.
Thus, ORD(t) becomes {1,1,0,0,1,1} after applying a logical-OR operation to RDL, and WL(t) becomes {3,5,0,0,6,1} by summing up the corresponding sensed raw data in RDL. Here, in the case of S = 1, we cannot estimate human movements actually. However, we can detect if humans exist in the detection area of the IR sensor. It is often called "proximity" position detection. To simplify the discussion, we treat S = 1 case as a human-trackable case in this paper.

Figure 2.
Example of sensed raw data processing.
Next, using ORD(t) and WL(t), the fired sensors during time period [t, t-WS1 + 1] are grouped using a clustering method. Finally, the weighted centers of the clusters are calculated, and they are used in the path estimation step as the candidates for the next human positions. The clustering algorithm is shown in Figure 3. The algorithm is based on the Ward method [26]; however, the termination condition of our algorithm is different from that of the original method. While the original Ward method tries to merge the nearest pair of clusters until all initial clusters are merged into one cluster, our algorithm stops the cluster merging if Equation (3), i.e., the termination condition, is no longer satisfied. The geometrical distance between the centers of any two clusters in the finally obtained cluster set will be longer than the radius of the detection range of IR sensor r.
Here, we assume that a person exists in each obtained cluster, and his/her position is the weighted center of the cluster. The weighted center coordinate cci(t) = (xi(t), yi(t)) of cluster ci(t) at time t is calculated using Equation (4).
Note that each element can belong to more than one cluster in our new clustering method, contrary to that used in our previous algorithm [24]. Thus, we can manage the case that multiple humans exist in the detection range of the same IR sensor, which cannot be managed by our previous algorithm introduced in [24 is the coordinate of element k , and | | c is the number of elements in cluster c: where "s" is a sensor ID. The weighted centers of all clusters are set to CC(t). They become candidates for human positions at time t. Figure 4 shows an example. In the example, cluster c1(t) consists of sensors s = 1, s = 2, and s = 3. If the position coordinates of the sensors are (10, 10), (9,5), and (14,8), and their weights are w1(t) = 4, w2(t) = 1, and w3(t) = 2, the center coordinate cc1(t) will be (11.0, 8.7), i.e., x1(t) = (4 × 10 + 1 × 9 + 2 × 14)/(4 + 1 + 2) and y1(t) = (4 × 10 + 1 × 5 + 2 × 8)/(4 + 1 + 2). Furthermore, sensor s = 2 belongs to another cluster c3(t). Thus, the coordinate of s = 2 is also used to calculate the weighted center cc3(t) of cluster c3(t). If a cluster contains only one sensor, like cluster c2(t) in Figure 4, the center position coordinate is simply the position of the sensor. After all weighted center coordinates of the clusters are calculated and set to CC(t), the path estimation step will be initiated.

Path Estimation Step
The procedure of this step is described in detail in Figures 5 and 6. In the main procedure, PathEstimation(), if TDL is empty, NewTarget() is called to create new humans. Their starting points  (1) and (2) are simply the weighted center coordinate of CC, respectively. Here, TDL represents humans currently existing in the room. If some humans exist, i.e., TDL is not empty, TrackTarget() in Figure 6 is invoked to obtain the next position coordinate of each human. TrackTarget(), which has been changed from our previous one introduced in [24], contains two sub-processes: prediction and modification. In the prediction sub-process, a predicted coordinate, td.PC, is provided for each human as the next position coordinate at time t (see "Block 2" in Figure 6). First, a movement angle α of each human at time t is calculated using a movement vector MV(td.ID). Here, the movement vector represents the expected direction of the corresponding human, and is calculated using the most recently estimated WS2 number of the position coordinates of the human, i.e., mvt(td.ID). Here, equal values for mcWS2(td.ID) and mc1(td.ID) indicate that the target human did not move. Therefore, the predicted coordinate td.PC is simply treated as mc1(td.ID). Otherwise, the prediction is calculated. Figure 7 shows an example. First, an average of movement distance mdAVE(td.ID) is calculated using mvt(td.ID), defined as: where x'i = mci(td.ID).x and y'i = mci(td.ID).y. Using the movement angle α and mdAVE(td.ID), the predicted coordinate td.PC is calculated (see also procedure for TrackTarget() shown in Figure 6).   Next, as the modification sub-process, all pairs of the weighted coordinates in CC and target humans in TDL are examined, and the best pair is selected in consideration of distance d, angle θ, and time to live (TTL) of "td." Here, "d" is the distance between td.PC and cc. "θ" is the relative angle between "mc1(td.ID) → td.PC" and "mc1(td.ID) → cc." A pair that minimizes "d" and "θ" and maximizes TTL at the same time is finally selected as the best pair.
Actual modification is, then, applied for the best pair found in the prediction sub-process, and the next movement position coordinate is obtained as the center of td.PC and cc of the best pair (see also "Block 3" in Figure 6). The obtained next position is appended to the td.Path of the corresponding target human. The elements representing the selected pair of the weighted center coordinate cc and target td are removed from CC and TDL. The aforementioned best-pair selection process is repeated until no more selected pairs exist. In addition, if the next position of a human can be successfully found, the TTL value of the human is incremented; otherwise, it is decremented. Here, the initial value of TTL is "1," and its maximum value is TTLMAX. The TTL value never exceeds TTLMAX. If the TTL of a human becomes less than 0, the human is removed from the TDL as he/she has left the room. In addition, if an unselected coordinate is contained in CC, NewTarget() is invoked and a human is appended to TDL as he/she enters the room. Due to this mechanism, even if the number of humans is not estimated correctly using the clustering algorithm in the location estimation step shown in Figures 1 and 3, and some humans have accidentally been generated by NewTarget() at this point, they will be removed within a short period. Figure 8 shows an example of the Path Estimation step. Here, TDL = {td1, td2} and CC = {cc1(t), cc2(t), cc3(t)}. First, a pair of td1.pc and cc1(t) is selected. The distance between them is smaller than that between others, and the angle is less than 180° and smaller than that between others. Thus, the center coordinate of td1.PC and cc1(t) is calculated, and it is appended to td1.Path. In addition, td1.TTL is incremented. In the next step, because TDL = {td2} and CC = {cc2(t), cc3(t)}, a pair of td2.PC and cc2(t) is selected, and the process described above is applied to the pair. Our path estimation algorithm shown in Figure 6 can be formulated using Equations (6) and (7). Assume that p is the position of a target human at time t, p = x y . Then, a prediction at time t is given by: where p is a predicted coordinate at time t; in this paper, it is the same as td.PC: p = x y . In addition, B is the average of the movement distance mdAVE(td.ID) of the target human, and u is the x-y component form of the predicted movement distance at time t. u = sin cos . Here, α is the angle of the movement vector MV(td.ID) at time t. Furthermore, the modified position at time t is given by: where z is a measurement of the target human at time t and is the same as the center position coordinate cc(t), which is obtained by the best-pair selection. z = x y . Equation (6) is similar in form to the Kalman filter, which is used in related work [19,20]. However, as compared to the original Kalman filter, in our proposed algorithm, p and z are calculated as mentioned previously, not using the method for obtaining the Kalman gain. Since the system environment is drastically changed by the numbers of IR sensors and humans in a room and the detection range and distribution density of IR sensors, we cannot apply the calculation method to obtain the Kalman gain.

Evaluation Environment
To evaluate our algorithm, we used generated simulation data. The details are as follows. First, the density of IR sensors "D" is defined by: where "S" is the number of IR sensors; "r" is the radius of the IR sensor detection range, r = 2.0 m; and "A" is the entire area of the room monitored, A = 100 m 2 (10 m × 10 m). We assume that the entire area A should be covered with the minimum number of IR sensors. Thus, if we fix the density D as 2.0, 3.0, 4.0, or 5.0, the number of IR sensors S that should be deployed is calculated using Equation (8). The corresponding value of S is actually 16, 24, 32, or 40 as shown in Table 1. Here, "S" sensors should be randomly deployed in the area.  The number of humans "H" is changed from one to four, and the human movement data are created for each case. The human movement scenario is as follows. After a random waiting time from 0 to 30 s, each human should enter an area through a door and walk around randomly. The walking scenarios are as follows. The human decides a destination position and walks toward it using the shortest path. Then, the human either stops at the destination position, or walks to a new random destination. Here, the human does not return quickly. After 60 s, the human should exit through the door. The walking speed of each human is randomly changed from 1.25 m/s to 1.75 m/s. For each "H," 50 different variations of data were generated. In total, 800 (4 different H × 4 different S × 50 different variations) different input data were generated and used for the evaluation of the proposed algorithm. Here, the parameters of our algorithm are WS1 = 4, WS2 = 24, and TTLMAX = 18. According to our preliminary estimations, these combinations of values produce good path estimations.

Evaluation Results and Remarks
Our method estimates the number of humans and their movement paths simultaneously. Thus, it is difficult to evaluate the method using tracking errors like in other human tracking methods, which assume that the number of humans is fixed and known [20,27]. To evaluate our method, we introduced four metrics: success estimation rate, averaged error, averaged tracking rate, and success rate of the number of humans. The details are as follows: Tables 2-5 show the evaluation results. For comparison, we also evaluated the previously proposed algorithm [24]. In the tables, the upper values represent the results of the proposed algorithm and the lower ones represent the results of the previous one.
( 1) Success estimation rate Table 2 represents the percentage of 50 input patterns provided for each different case that were estimated correctly by our algorithm. For H = 1, both algorithms estimated all human movement paths correctly for all different densities. For H = 2, the success estimation rate of the proposed and previous algorithms were 60% and 57% on average, respectively. In addition, for H = 3 and H = 4, both the estimation rates were decreased to less than 34%. If the number of humans increase, their movement paths are often overlapped, which makes it difficult to estimate individual human movement. Thus, the success estimation rate decreases with the increase in the number of humans in a room, regardless of the path estimation methods being compared.  (2) Averaged error Table 3 shows the averaged error of the estimated human movement paths compared to the true human movement paths among the successfully estimated data. In our algorithm, a position coordinate of the IR sensor was basically used. Therefore, an error in the detection range of the IR sensor is contained during the initial stage of the path estimation. Even so, for H = 1 and H = 2, the proposed method can trace the human movement path accurately for less than 0.60 m and 1.42 m on average, respectively. However, although our previous algorithm can trace one human for less than 0.68 m, the averaged error is more than 2.08 m when H = 2. In other cases, it could not trace the human movement paths correctly. (3) Averaged tracking rate In Table 4, the averaged tracking rate for each target human is described. The tracking rate is defined as the ratio of the length of the correctly estimated route to the total length of the actual human movement path. If the rate is high, it indicates that the human path is estimated well for a long time. For H = 1, the human was tracked with an accuracy rate of more than 94% in both algorithms. In other cases, if the density was high, the tracking rate was also high. This implies that many sensors are needed to estimate multiple humans.  (4) Success rate of the number of humans Table 5 shows what percentage of the number of humans correctly estimated in the room at each time. It is directly related to our clustering algorithm. Using the proposed algorithm, when H = 1 and H = 2, the success rates of H were more than 97% and 76%, respectively. However, for H = 3 and H = 4, these rates were from 45% to 60%. If the number of humans increases, the probability that some of them are walking close together increases. Therefore, the performance of the clustering algorithm decreases. Compared to our previous algorithm, the newly proposed clustering algorithm does not perform well. According to these results, if the number of humans in the room increased, the tracking accuracy of our method decreased. This is because our clustering algorithm does not work well if some humans walk closely together. In addition, if relatively many humans exist, e.g., when H = 3 and H = 4 in our experiments, almost all the IR sensors are fired simultaneously. In such cases, no algorithm can estimate the number of humans correctly. However, as shown for H = 2, our algorithm can estimate the human movement paths well with a 1.42 m averaged error. In addition, our human-location-estimation algorithm based on a clustering method performs with a success rate of greater than 76% even for the complicated case when two persons enter and exit the room at different times. Compared to our previous algorithm, the clustering method itself in the proposed algorithm does not perform well, because the new clustering method allows a position candidate to belong to more than one cluster. However, owing to this fact, the averaged accuracy of the final human tracking is improved. As shown in Table 3, the averaged errors of the estimated human movement paths are improved two-fold for H = 2. Figure 9 illustrates an example of the estimated human movement paths for H = 2. The solid line is the true path of the first human, and the dotted line is the true path of the second human. They are indicated by "TruePath1" and "TruePath2" in Figure 9. Two different kinds of markers represent the estimated locations of two persons, respectively. They are "EstPath1" and "EstPath2" in the figure. The start and goal locations of both persons are the same, which is indicated as "Door" in the figure. As shown in Figure 9, our proposed algorithm can trace human movement well. This is the best result among 50 different input data generated with the same condition, D = 5.0. That is, 40 sensors were deployed in the field. The numerical results for the case shown in Figure 9 are as follows: the averaged error of the estimated paths, corresponding to Table 3, is 0.56 m. The averaged tracking error, corresponding to Table 4, is 99.89%. Success detection rate of number of humans, corresponding to Table 5, is 86.85%.

Conclusions
We have proposed an algorithm that can track the human movement paths using only the binary sensing data obtained from infrared sensors attached to the ceiling. The human positions are estimated at each timepoint based on a clustering method. Thus, the proposed algorithm can track multiple humans even if the number of humans in the room is changed dynamically, which was difficult to realize using methods proposed in related studies. According to simulation-based evaluations, our algorithm can trace real human movement paths with a 1.32 m error on average if two humans are in the room. In future studies, we will evaluate our algorithm using real sensed data obtained from a real infrared sensor system.


Exploring partitioning methods for multicast in 3D bufferless Network on Chip

In this paper, we first propose two multicast partitioning methods named TBP (two block partitioning) and LBP (layer block partitioning), each of which has a different level of efficiency in 3D bufferless network on chip. In addition, we present a Recursive Partitioning (RP) method in which the network is recursively partitioned until all partitions contain comparable number of multicast destination nodes. By this method, the multicast destination nodes are distributed evenly and the network latency is significantly decreased. Simulation results illustrate that the RP scheme achieves 39%, 45% and 41% less latency on average than that of the TBP scheme and 8%, 17% and 12% less latency on average than that of the LBP scheme under three synthetic traffic patterns respectively.

Introduction
As technology scaling with Moore's law, a single die can integrate several hundreds of cores. When the number of cores increases, interconnect becomes more and more complicated in System-on-Chip design. Networks on chip (NoC) were shown to be feasible and easy to scale for supporting a large number of cores and memory [1]. NoC can provide high throughput and low latency for one-to-one (unicast) traffic at present. In order to run many kinds of parallel applications on Multi-Processors System-on-Chip (MPSoC), collective communications must be supported. In an MPSoC system, the cache-coherent shared memory protocols (such as directory-based or token-based) require one-to-many or broadcast communications to obtain shared data or invalidate shared data on different cache blocks [2]. It has been proved that only 2-5% multicast traffic in the total network traffic will have a serious impact on system performance [3]. As a result, supporting efficient multicast in NoC can make the performance of the MPSoC system better.
Furthermore, with the increasing number of cores, the 2D plane is not an efficient way due to its long wire interconnections [4]. The emerging of the 3D integrated circuit stacks several dies to reduce the long wire delay, which results in lower power consumption and higher performance [5,6]. The 3D integrated circuit is also an attractive way for Network on chip. The 3D NoC architecture is widely studied in the network topology [7], router architecture [8], and routing algorithms [9].
Previous researches on NoC focus on the wormhole/virtual-channel router which contains buffers in each input port to buffer the packets. Though buffers can improve the bandwidth efficiency of the network, they also consume much energy and chip area [10]. Bufferless NoC has become a potential alternative to the traditional wormhole based NoC lately [10,11,12]. The bufferless router makes use of deflection routing to eliminate the buffers, which has two benefits: simplicity in design and reduce hardware cost. However, in bufferless NoC, packets cannot wait in a router and deflection routing makes the routing path unpredictable. Therefore it is impossible to employ existing multicast algorithms used in networks with buffers.
A method related to our work is discussed in [13] and [14] which propose partitioning methods for the 3D mesh network with buffers. However, it is not suitable for 3D bufferless network due to the fact that unpredictable routing path exists in bufferless network. To the best of our knowledge, there is no work addressing multicasting in 3D bufferless NoC systems. In this paper, we first propose two multicast partitioning methods named TBP (two block partitioning) and LBP (layer block partitioning), each of which has a different level of efficiency in 3D bufferless network on chip. In addition, we present a Recursive Partitioning (RP) method in which the network is recursively partitioned until all partitions contain comparable number of multicast destination nodes. By this method, the multicast destination nodes are distributed evenly and the network latency is significantly decreased. Simulation results illustrate that the RP scheme achieves 39%, 45% and 41% less latency on average than that of the TBP scheme and 8%, 17% and 12% less latency on average than that of the LBP scheme under three synthetic traffic patterns respectively.

3D bufferless router architecture
The baseline 3D bufferless router is a 1-cycle router with a 3-stage permutation network [8]. The structure of the baseline 3D bufferless router is show in Fig. 1. There are 7 input/output ports in each router. Each input port has only one input register, therefore the packet is not buffered in the router. It utilizes deflection routing to route packets. In the case of two or more flits contending for the same output of the permutation cell, the permutation cell selects the flit with the higher priority to the shortest output port and misroutes the other flit to the other output. The priority is decided according to the number of hops the packet has been routed in the network which can limit the number of misroutings to avoid livelock. The router can handle at most six packets at the same cycle and make output port allocation for each packet from high priority to low.

Partitioning methods for 3D bufferless mesh architecture
In this section, we propose three different partitioning methods for 3D bufferless mesh architecture.

Two-block partitioning (TBP) method
As the number of multicast destination increases, the long path of baseline 3D bufferless deflection routing without packets replication will lead to large multicast latency. To solve this problem, we propose TBP scheme with adaptive packets replication at the source node. Fig. 2 shows an example of this partitioning policy and the portions of each partition depending on the source node position. It is assumed that the source node is node 25 while the destination node set is D ¼ f9; 16; 17; 24; 39; 44g. The multicast based on TBP method is shown in Fig. 2

(a).
First the multicast packet in Node 25 is copied into two and sent to the two different partitions separately, if the source router has more than two empty ports. Node 24 and 39 are chosen as the first best candidate since they have the minimum Manhattan distance to the source node 25. After the packet is sent to node 24 and 39, node 16 and 9 are chosen as the second best candidate. Without contention, the minimum multicast latency is equal to 9 hops. The path shown in Fig. 2(a) is not the only one path since the packet may be deflected due to contention. As illustrated in Fig. 2(a), the source node is located at the middle layer, two partitions cover comparable number of nodes. However in Fig. 2(b), the source is located at the first layer, one partition contains much more nodes than the other. Now suppose that the multicast message m ¼ ð6; f1; 16; 17; 24; 39; 44gÞ is generated by the source node 6. The destination nodes are split into two groups according to the labels: G H ¼ f16; 17; 24; 39; 44g and G L ¼ f1g. The packet created for G H uses the manhattan path as follows: f6; 25; 24; 39; 40; 47; 16; 17; 46; 45; 44g where 10 hops are needed to reach the last destination. The packet path for the G L is: f6; 1g only one hop is required for delivering the packet to the destination. This is the drawback of TBP method.

Layer block partitioning (LBP) method
In this method, the network is partitioned into layer 0; 1; . . . ; n subnetworks. Destination are sorted based on each layer. The nodes with the same z value will be put into the same group. The scheme has several advantages over the TBP scheme as it achieves a high degree of parallelism and avoids the creation of long multicast paths. For the LBP scheme, when a multicast packet is being injected into the network from the local port, if free output ports are plane ports, the packet will not be replicated and routed the same way as the path-based 3D bufferless multicast scheme. If the number of free output ports is more than one and is the up or down port and destinations fall into different layer corresponding to the free output ports, the destinations will be grouped according to the layer partition and the multicast packet will be replicated and routed through different output ports. The pseudo code of LBP scheme is shown in Table I.
Free_port is a 6-bit vector, each bit of which indicates the corresponding output port is free or not. Output_prio sorts the priority of the output port from the highest priority to the lowest. The router always chooses a free output port with the highest priority to replicate and delivery a copy of the multicast packet. The function replicate_addr_gen acquires the address of the replicated packet. After replication, the original destination address will be updated and the number of packet copies will increase by 1. Fig. 3 shows an example of the multicast routing algorithm based on the LBP method. For the multicast message m ¼ ð25; f9; 16; 17; 18; 24; 39; 46gÞ, three groups are formed: layer0 ¼ f9g, layer1 ¼ f16; 17; 18; 24g, layer2 ¼ f39; 46g. Without contention, the multicast packet requires 5 link traversals in total. The path shown in Fig. 3(a) is not the only one path since the packet may be deflected due to contention. The unbalanced partitions are shown in Fig. 3(b). The multicast message m ¼ ð6; f1; 9; 16; 17; 18; 24; 39; 46gÞ requires 6 hops in total to reach the final destination without contention.

Recursive partitioning (RP) method
The objective of the recursive partitioning method is to optimize the number of destination nodes that can be included in a partition. In this method, the network is recursively partitioned until each partition includes less than n destination nodes. In the worst case, the network is partitioned into c layer partitions like the LBP method in an a Â b Â c network. So, we have considered the value n as the maximum number of destination nodes in a partition of the LBP method. The pseudo code of RP scheme is shown in Table II. Num_dst indicates whether a partition of the LBP method includes more than n destination nodes. The number of n will change if the number of destinations of the multicast packets is changed or the size of the network is changed. When free_port and num_dst are equal to 1, the router will choose a free output port to replicate destination address and delivery a copy of the multicast packet. The function replicate_addr_gen_1 acquires the address of the replicated packet in a partition of the LBP method. After replication, the original destinaiton address will be updated and the number of packet copies will increase by 1.
An example of the multicast routing algorithm based on the RP method is illustrated in Fig. 4, where a multicast message is generated at the source node 25. The steps of the RP algorithm can be expressed as follow: 1. The value n is set to 4 in a 4 Â 4 Â 3 network; 2. The network is divided into three partitions using LBP method. The Fig. 4(a) shows the partitions when the source node is 25; 3. If the number of nodes in a partition exceeds the predefined value n, the function replicate_addr_gen_1 acquires the address of the replicated packet in a partition of the LBP method. the partition is divided into new partitions according to the location of the destination nodes. A region partition policy is proposed to divide the layer into at most 4 regions (P1,P2,P3,P4) according to the position of the source node, as shown in Fig. 4(a) layer1. If there is only one free output port, the packet will not be replicated and will be routed the same way as the LBP scheme. If the number of free output ports is more than one and destinations fall into different regions corresponding to the free output ports, the RP scheme will work. Fig. 4(a) shows the example of recursive partitioning. There are 5 destination nodes in layer1 which is greater than n = 4. So, the multicast packet will be replicated according to the region partition in layer1. Without contention, the multicast path requires 5 hops. However, the multicast path will require 12 hops when it using the LBP method.
The same partitioning technique is applied to the low layer subnetwork. Fig. 4(b) shows another example of the RP method where the source node is 6. In this example, the multicast path requires 6 hops as the maximum latency.

Experimental evaluation
In this section, we evaluate the performance of the three proposed multicasting partitioning mechanism for the 3D mesh interconnection network with a cycleaccurate NoC simulator which is developed in VHDL.

Methodology
The experiments are performed on a 4 Â 4 Â 4 3D mesh network. The routing process is partitioned into two steps. First of all, the routing computation in company with the input and output priority sorting are executed in parallel. Then the output port allocator allocates each packet based on the results of the first step. For the LBP and RP routers with adaptive packets replication, multicast packets are replicated at the moment of output allocation when at least one free output port left for replication. The three multicast routers are more complicated than the regular unicast bufferless router as a result of increasing route computation for multicast packets and packets replication. For the purpose of simplify the design, we use a 1cycle router [8] for three kinds of multicast routers (TBP, LBP, RP routers) which is shown in Fig. 1. To perform the simulations, a combination of unicast and multicast traffic is used. For unicast traffic, we use three synthetic traffic patterns: uniform random, bit complement, transpose. For multicast packet, the destination positions are uniformly distributed.
We measure the average packet latency T which is calculated by the equation and measured in cycle, where T src is the time that a packet startup latency and T net is the network delivery time. We also evaluate the link utilization, which is the average percentage of used links in total links. Link utilization can reflect the power consumption. The more used links, the higher power consumption is consumed.

Latency analysis
We measure the performance of the three multicast schemes in this subsection. Fig. 5(a)-(c) illustrate the average packet latency of the three multicast schemes. The multicast traffic is 10% of the total traffic and the number of destinations is 8. It can be observed that the RP scheme achieves the smaller average latency than the other two multicast schemes. However, the network with the RP scheme reaches saturation earlier than the other two schemes. As the network is not saturated, the average latency of the LBP scheme is 32%, 40% and 35% less than that of TBP scheme under three synthetic traffic patterns respectively. The RP scheme achieves 39%, 45% and 41% less latency on average than that of the TBP scheme and 8%, 17% and 12% less latency on average than that of the LBP scheme under three synthetic traffic patterns respectively. The results are consistent with our expectations. The RP scheme outperforms the LBP scheme and TBP scheme at low or medium traffic loads. The RP scheme performs similar as the LBP scheme when the network reaches saturation point. The fact is that they can replicate multicast packets at both source and intermediate nodes which increases the network loads, so the throughput of the network is a little bit lower than TBP scheme. However, under real application workloads the network is never fully loaded, so the RP scheme is a better alternative for bufferless NOC than the other two schemes. We synthesize three routers using Synopsys Design Compiler with 65 nm technology at the typical condition (25°C and 1.0 v). The critical path delay of the TBP multicast routers is 2 ns while the LBP and RP methods are at most 8% higher than this value, thus all the above conclusions still hold. Fig. 6 shows the link utilization with the number of multicast destinations increasing from 4 to 16 under random traffic with 10% multicast traffic. In the case of the network not saturated, compared with TBP scheme, the RP scheme can save 12%, 16% and 23% of link utilization on average for different number of multicast destinations respectively. Compared with LBP scheme, the RP scheme can save 3%, 5% and 10% of link utilization. It can also be concluded that as the number of destinations increases, the RP scheme can save more power consumption than that of the other two multicast schemes.

Scalability
We also evaluate the performance of the three multicast schemes with different portion of multicast traffic and different number of multicast destinations. Fig. 7(a) shows the performance variation of the three multicast schemes with different  portion of multicast traffic from 10% to 40%. The number of multicast destinations is 8. It can be seen from the figure, the RP scheme achieves much less average latency than the TBP and LBP schemes with the increasing portion of multicast traffic. Fig. 7(b) reveals the average latency of the three multicast schemes with the increasing number of multicast destinations. The portion of multicast destinations is 10%. As the number of multicast destinations increases from 4 to 32, the average latency of the RP scheme does not vary so much, while the average latency of the LBP scheme increases by up to 1 times and the average latency of the TBP scheme increases by up to 3 times.

Conclusion
This paper proposes a novel idea of balanced partitioning which leads the 3D mesh architecture to be partitioned effectively. We first presented TBP and LBP partitioning methods for 3D bufferless mesh architecture, each of which has a different level of efficiency. In order to achieve higher performance, the recursive partitioning (RP) method was introduced. This method partitions the network recursively until all partitions contain considerable number of nodes. Our simulations with different traffic profiles showed that the RP method reduces the average packet latency and provides higher degree of parallelism compared with the other two methods.