Hoang Thanh Lam

Wenjie Pei

Adriana Prado

Baptiste Jeudy

Élisa Fromont

Mining Top-K Largest Tiles in a Data Stream

Large tiles in a database are itemsets with the largest area which is defined as the itemset frequency in the database multiplied by its size. Mining these large tiles is an important pattern mining problem since tiles with a large area describe a large part of the database. In this paper, we introduce the problem of mining top-k largest tiles in a data stream under the sliding window model. We propose a candidatebased approach which summarizes the data stream and produces the top-k largest tiles efficiently for moderate window size. We also propose an approximation algorithm with theoretical bounds on the error rate to cope with large size windows. In the experiments with two real-life datasets, the approximation algorithm is up to hundred times faster than the candidate-based solution and the baseline algorithms based on the state-of-the-art solutions. We also investigate an application of large tile mining in computer vision and in emerging search topics monitoring.

Introduction

Mining frequent patterns is an important research topic in data mining. However, instead of focusing on exhaustive search to find all possible frequent patterns, many works are now focusing on designing methods that are not only efficient in the context of very big data but also limit, e.g. with constraints or with new interestingness measures, the number of patterns output by these algorithms.

Area is a measure of pattern interestingness defined as a pattern's frequency in the database multiplied by its size. It has been shown that in some applications such as in role mining [START_REF] Lu | Optimal boolean matrix decomposition: Application to role engineering[END_REF][START_REF] Vaidya | The role mining problem: A formal perspective[END_REF] where the idea is, given a set of users and a set of permissions, to find a minimum set of roles such that all users will be assigned a role for which some permissions will be granted, mining roles is equivalent to a variant of mining itemsets with large area in a database.

Recent applications produce a large variety of transactional data streams, such as text stream from twitter 4 or video stream in which video frames can be converted into transactions [START_REF] Fernando | Effective use of frequent itemset mining for image classification[END_REF]. In the context of data streams, data usually arrive continuously with high speed, hence requiring efficient mining techniques for summarizing the data stream and keeping track of important patterns.

In this paper we tackle the problem of mining the top-k largest tiles, i.e. the k closed itemsets with the largest area in a stream of itemsets. The problem of mining the largest tile in a database is well-known to be NP-hard and inapproximable [START_REF] Geerts | Tiling databases[END_REF]. Therefore, the straightforward approach that recalculates the set of top-k largest tiles from scratch every time a sliding window is updated is not efficient. To deal with this situation, we first introduce in Section 4 a candidate-based approach which incrementally summarizes the stream by keeping the itemsets that can be the top-k largest tiles in some future windows. For each candidate, an upper-bound and a lower-bound of the area in any future window are kept. These bounds are used to prune the candidates when they cannot be the top-k largest tiles in any future window. In doing so, the candidate-based algorithm is more efficient than the straightforward approach because updating is cheaper than recalculating the top-k tiles from scratch.

However, when the widow size is large, the candidate-based algorithm is inefficient because the summary grows very quickly. Therefore, we introduce an approximation algorithm with theoretical bounds on the error rate. In the experiments with two real-life datasets presented in Section 6, the approximation algorithm is two to three orders of magnitude faster and more accurate than the candidate-based solution and the baseline algorithms created based on the stateof-the-art solutions for the top-k largest tiles mining. We also discuss potential applications of mining large tiles for object tracking in video and emerging search topics monitoring problems.

Related Works

Recent works [START_REF] Pasquier | Discovering frequent closed itemsets for association rules[END_REF][START_REF] Van Leeuwen | Diverse subgroup set discovery[END_REF][START_REF] Vreeken | Krimp: mining itemsets that compress[END_REF] propose approaches to solve the redundancy and trivial patterns issues in frequent pattern mining. For instance, the first two works focuse on finding a relevant or concise representation of sets of frequent patterns. Meanwhile, the last work solves the aforementioned issues by proposing a MDLbased approach that finds patterns compressing the database well. In all these cases as well as in ours, the purpose is to limit the output of the algorithms to a set of useful patterns.

The problem of mining large tiles has already been tackled in [START_REF] Lu | Optimal boolean matrix decomposition: Application to role engineering[END_REF][START_REF] Vaidya | The role mining problem: A formal perspective[END_REF][START_REF] Remmerie | Unraveling tobacco by-2 protein complexes with {BN} page/lcms/ms and clustering methods[END_REF]. The authors of the two first papers show that the problem of finding roles is equivalent to variants of the tiling database problem. The last paper shows how tiles can be used to understand complex proteins by identifying their subunits. The authors of [START_REF] Cerf | Closed and noise-tolerant patterns in n-ary relations[END_REF] also showed that the dense rectangles in a binary matrix seem to correspond to interesting concepts. Therefore, mining large tiles may help to identify those interesting concepts from the databases. In [START_REF] Tatti | Discovering descriptive tile trees by mining optimal geometric subtiles[END_REF] the authors investigate how to output the set of tiles in a tree representation which is easily interpretable by the users. Tiles are also used to identify and characterize anomalies in a database in [START_REF] Smets | The odd one out: Identifying and characterising anomalies[END_REF].

In many applications, data arrives in a streaming fashion with high speed [START_REF] Aggarwal | Data Streams -Models and Algorithms[END_REF]. Besides the popularity of data streams in many applications, the temporal aspect of the patterns such as the evolution of patterns overtime [START_REF] Aggarwal | Data Streams -Models and Algorithms[END_REF] provides

a b c d e f g h k 𝑇 1 0 0 0 1 1 0 1 1 0 𝑇 2 1 1 1 1 1 0 1 1 0 𝑇 3 1 1 1 1 1 1 1 1 1 𝑇 4 1 1 1 0 0 0 1 1 1 𝑇 5 1 1 1 0 0 0 1 1 1 𝑇 6 1 1 1 0 1 1 0 0 0 𝑇 7 0 1 1 0 0 0 1 1 1 𝑇 8 0 0 0 0 1 1 1 1 1
Fig. 1. An example of large tiles in a window with w = 8 transactions.

useful insights about the data for the users. Despite the importance of data stream paradigm, there is no work yet addressing the problem of mining large tiles in a data stream. The algorithms introduced in this paper are, to the best of our knowledge, the first to solve the problem of mining the top-k largest tiles in a stream of itemsets under the sliding windows model.

Problem Definition

Let Σ be an alphabet of items, a transaction T is an itemset

T ⊆ Σ. Let S = T 1 T 2 • • •
T n be a stream of transactions where each transaction T t is associated with a timestamp t indicating the order of transaction generation.

In this work we consider the sliding window model in which a window of the w most recent transactions is monitored. A sliding window of size w at time point t denoted by W t is a sequence of w consecutive transactions in the data stream, i.e.

W t = T t-w+1 • • • T t-1 T t .
For any given itemset I ⊆ Σ the frequency of I in a sliding window W t , denoted by f t (I), is defined as the number of transactions in W t that are the supersets of I. Let |I| be the cardinality of the set I, the area of I in the window W t denoted by A t (I) is defined as A t (I) = f t (I) * |I|.

An itemset I is closed in a window W t if there is no superset of I with the same frequency as I in the window W t . Such itemsets are usually called large tiles in the literature [START_REF] Geerts | Tiling databases[END_REF]. Itemsets that are not closed correspond to sub-tiles of large tiles and are thus not interesting for our problem. From now on we use the term tiles to refer to closed itemsets.

Example 1 (Large tiles). Figure 1 shows a window with 8 transactions represented as rows of a binary matrix. In each row, an element is equal to 1 if the corresponding item belongs to the transaction. In this figure, three large tiles abc, degh and ghk with respective area 15, 12 and 15 are highlighted. Large tiles are the maximal sub-matrices containing only 1.

The problem of stream tiling can be formulated as : Definition 1 (Stream Tiling). Given a data stream of itemset transactions S, a parameter k and a window size w, the stream tiling problem consists in computing the k tiles with the largest area in every window of size w.

Algorithms

It was proven that the problem of mining the largest tile in a database is NP-Complete and is even inapproximable [START_REF] Geerts | Tiling databases[END_REF]. Therefore, recalculation of the topk largest tiles from scratch every time the window is updated is very timedemanding. In this section we discuss efficient solutions for the given problem under the streaming context.

Candidate-Based Algorithm

We first discuss an exact algorithm named cTile which maintains a summary of the sliding window containing candidate itemsets that potentially can become top-k largest tiles in any future window. This summary is designed such that it can avoid expensive recalculation of the set of largest tiles from scratch when the window is sliding.

The general idea of the algorithm is as follows: at every transaction T i we keep a candidate list C i of all closed subsets of the given transaction which can become a top-k tile in a future sliding window. In order to identify these closed itemsets, we keep a lower-bound and an upper-bound on the area of these itemsets in any future sliding window which contains T i . These bounds will be used to infer which itemsets can be top-k tiles and which sets for sure cannot be top-k tiles in any future window (thus can be removed from the summary).

𝑇 1 𝑇 2 … … … … 𝑇 𝑖 … … … 𝑇 𝑤-1 𝑇 𝑤 𝑊 𝑤 𝑊 𝑤+𝑖-1 𝐽 -= min 𝐴 𝑤+𝑖-1 (𝐽) 𝐽 + = max 𝐴 𝑤+𝑖-1 (𝐽) 𝑖 -1 transactions Given a window W w = T 1 • • • T w-1
T w , for every closed itemset J in the candidate list C i (1 ≤ i ≤ w), the lower-bound and the upper-bound of the area of this itemset in any future sliding window containing T i are denoted J -and J + . The lower-bound J -at time point w is calculated as the area of J in the transactions T i , T i+1 , • • • , T w and the upper-bound J + is calculated as

J + = J -+ (i -1) × |J| (1)
Proposition 1. J -and J + are the correct lower-bound and upper-bound on the area of J in any future window W t of size w, w ≤ t, containing the transaction T i . Algorithm 1 incrementally maintains the candidate lists C i in the summary. When a new transaction is appended to the end of the window, the bounds J + and J -are updated accordingly (line 10-13). Then the new transaction is intersected with every element in the summary and the intersections are added to the corresponding candidate list (line 8-9). If there exists another candidate B ∈ C j kept at a younger transaction T j (j > i), such that B -is ranked k th in all transactions younger than T i (including itself) and meanwhile B -> J + , then J will never be among top-k tiles in any future sliding window, hence J is removed from the summary. When a transaction is expired, it is removed from the window along with the candidate list stored at the given transaction. After that all the lower-bounds and the upper-bounds are updated accordingly.

Example 3 (cTile)

Assume that we want to get top-3 largest tiles in the sliding window, i.e. k = 3. Since (b, 3, 3) in C 2 has the upper-bound equal to 3. It is removed from the summary because the itemset ranked at the third position in the candidates lists C ← C \ {Ct-w+1} 6:

Ct+1 ← {Tt+1} 7:

C ← C ∪ {Ct+1} 8:

for i = t → t -w + 2 do 9:
Ci ←-Ci intersect(Ci, Tt+1) 10:

for J ∈ Ci do 11:

J -= |{i ≤ l ≤ t|J ⊆ T l }| * |J| 12: J + = J -+ (w -t + i -1)|J| 13:
end for 14:

Pruning(Ci) 15:

end for 16: end for 17: Return C of younger transactions (including T 2 ) are (a, 4, 4), (ac, 4, 10), (ab, 4, 4) whose lower-bound is 4. The same pruning operation can be applied for (b, 2, 3) in C 3 .

Theorem 1 Given k, w and a stream of transactions, using the summary in Algorithm 1 we can answer the top-k largest tiles exactly.

Proof. Assume that there exists a top-k largest tile I in a window W t not recognized as a top-k largest tile (false negative). The first transaction of W t containing I is denoted by T i . False negative can only occur if at time t, I ∈ C i . Therefore either I is directly pruned from C i by pruning criteria or indirectly pruned because its closed supersets are pruned.

Since the bounds are exact the former case cannot happen. The latter case happens when I is not a closed itemset at the moment its closed supersets are pruned. In such case, the upper-bound of I is always less than the upper-bound of its closed supersets which were used to prune the supersets. Therefore, I is not a top-k tile in W t . Both cases lead to contradiction.

Approximation Algorithm

The size of the summary maintained by cTile grows quickly when the window size increases making it inefficient for monitoring large windows. The main reason is that each time a new transaction arrives, it has to be intersected with a large amount of candidates, which is a time-consuming operation. Therefore, in this section we discuss an approximation algorithm named aTile that approximates the set of largest tiles efficiently.

The main process of the approximation algorithm is almost the same as Algorithm 1. The only difference of these two algorithms lies in the method of candidate pruning. Instead of using an upper-bound on the area, aTile tries to approximate the future area of a tile. This estimate is used to prune the candidates in the same way as cTile does.

Let us denote the probability of observing an itemset J as a subset of a transaction in the data stream by µ J . Therefore, the expectation of the area of the itemset J in a window of size w is µ J .|J|.w.

Given a window W w = T 1 • • • T w-1 T w , assume that J is kept in the candidate list C i of transaction T i . The lower-bound J -is calculated as in the cT ile algorithm. Instead of using an upper bound J + on the area of J, we compute an estimate of the area of J denoted by J * .

Since µ J is unknown, we cannot directly define J * as µ J .|J|.w. The probability µ J is thus estimated in the sub-window T i • • • T w . However, if this sub-window is too small, the estimation of µ J is not accurate. Therefore, we introduced a threshold L ≤ w: if the sub-window is smaller than L, we do not use this estimation and fall back on using the upper bound J + for J * :

J * = J -+ (i -1) × |J| if w -i + 1 ≤ L (i.e., J * = J + ) J - (w-i+1) × w if w -i + 1 > L (2)
All the steps of the of the aTile algorithm are very similar to the cT ile algorithm except that it uses J * instead of J + for candidate pruning. In the experiment section we empirically show that the aTile algorithm is more efficient than the cT ile algorithm because its candidate set is more concise. An important property of the aT ile algorithm is that the error rate on the accuracy of the result is bounded if we assume that the transactions are i.i.d. in a window of size 2w.

Theoretical Analysis

In this section, we show theoretical bounds on the probability of errors induced by the aT ile algorithm. We mainly show that the error rate is extremely low when L is large enough. We consider two types of error event: False negative (FN): in a window W t * for some t * , there is a true top-k largest tile that is not present in the results returned by the aT ile algorithm and False positive (FP): in a window W t * for some t * , there is a non top-k largest tile that is present in the result returned by the aT ile algorithm.

False Negative Bound

Let us assume that I is a true top-k largest tile in the window W t * and I t * k is ranked at position k in the list of a true top-k largest tile of the window W t * . Let ≥ 0 be the difference between the average area of I and

I t * k in the window W t * , i.e. = A t * (I)-A t * (I t * k ) w
. The following lemma show the relationship between the probability of a false negative, L and :

Lemma 1 If the transactions are independent and identically distributed (i.i.d) in a window of size 2w, the probability that a random top-k tile is not reported, is bounded as follows:

P r(F N ) < 4w * e -L 2 2|I| 2
Proof. A false negative happens in the window W t * with respect to an itemset I if its area is underestimated in that window. This only happens if there exists at least one moment t < t * such that I is pruned from the candidate list C i (t -w < i ≤ t).

Given a time t, let W 0 be a window containing transactions T i , T i+1 , • • • , T t where |W 0 | = w o > L and let f 0 (I) be the frequency of I in the window W 0 . When t is given, the event "I is removed from C i " (denoted by R t ) only happens when the estimate of the upper-bound on the area of I is less than the area of

I t * k , i.e. f W0 (I)|I| w w0 < A t * (I t * k )
. Therefore, we have:

P r(R t ) < P r f 0 (I)|I| w w 0 ≤ A t * (I t * k ) < P r f 0 (I)|I| w w 0 < A t * (I) -w < P r f 0 (I) w 0 < f t * (I) w - |I| < P r f 0 (I) w 0 - f t * (I) w > |I| < P r f 0 (I) w 0 -µ I + f t * (I) w -µ I > |I| < P r f 0 (I) w 0 -µ I > 0 w 0 + P r f t * (I) w -µ I > 1 w < P r (|f 0 (I) -µ I w 0 | > 0 ) + P r (|f t * (I) -µ I w| > 1 )
Where 0 = w0 2|I| and 1 = w 2|I| . It is important to notice that w 0 µ I and wµ I are the expectation of the frequency of I in the window W 0 and the window W t * respectively. Since the transactions are independent to each other, according to the Hoeffding inequality [START_REF] Motwani | Randomized algorithms[END_REF] we have:

P r (|f 0 (I) -w 0 µ I | > 0 ) < 2e -2 2 0 w 0 < 2e -w0 2 2|I| 2 < 2e -L 2 2|I| 2 (3)
A similar inequality can be obtained for bounding the second term as follows:

P r (|f t * (I) -wµ I | > 1 ) < 2e -2 2 1 w < 2e -w 2 2|I| 2 < 2e -L 2 2|I| 2 (4)
Moreover, since F N happens only when there at least one moment t such that R t happens, therefore:

P r(F N ) < P r(∪ t * -w<t≤t * R t ) < t * -w<t≤t * P r(R t ) (5)
Inequalities 3, 4 and 5 prove the lemma.

A direct corollary of Lemma 1 is shown in the following theorem:

Theorem 2 If I is strictly more important than I t * k , i.e. the expectation of the area of I in a transaction is strictly greater than the expectation of the area of I t * k in a transaction and L = O(w) then: lim

L-→∞ P r(F N ) = 0
Proof. Let E(A(I)) be the expectation of the area of I in a transaction and E(A(I t * k )) be the expectation of the area of I t * k in a transaction. Since =

A t * (I)-A t * (I t * k ) w
we can imply that:

|I| = A t * (I) -A t * (I t * k ) w|I| E(A(I)) -E(A(I t * k )) |I| (6)
The last equation is the result of the law of large number when L goes to ∞.

From the last equation we can imply that lim

L-→∞ 4we -L 2
2|I| 2 = 0 from which the theorem is proved.

False Positive Bound

In this subsection we prove bound for false positive error which can be obtained in a similar way as the bound of false negative.

Let J be an itemset that is not a true top-k largest tile in the window W t * but returned by the aT ile algorithm as a false positive tile. An itemset ranked at position k in the list of a true top-k largest tile of the window W t * is denoted by I t * k . Let ≥ 0 be the difference between the area of J and I t * k in the window W t * , i.e. =

A t * (I t * k )-A t * (J) w . The following lemma show the relationship between the probability of a false positive and L, (proof is similar to Lemma 1):

Lemma 2 If the transactions are i.i.d. in a window of size 2w the probability of false positive, i.e. the event that a random non top-k tile is reported, is bounded as follows:

P r(F P ) < 4w * e -L 2 2|I| 2
A corollary of Lemma 2 is shown as follows (proof is similar to Theorem 2):

Theorem 3 If tile size is bounded and if J is strictly less important than I t * k , i.e. the expectation of the area of J is strictly less than the expectation of the area of I t * k in a transaction and L = O(w) then: lim L-→∞ P r(F P ) = 0 Theorem 3 and Theorem 2 show an interesting result that the probability of false positive and false negative decrease exponentially with L = O(w). It is important to notice that the condition L = O(w) can be replaced by a weaker assumption L = k log w for some constant values k. If k is large enough the bound is also closed to zero. Although the bounds are not tight, in experiments, we empirically show that false negative rate and false positive rate are negligible even when L is set to a small value.

Long Lasting Tiles

The two previous Theorems give probabilistic results. We can also show a deterministic one: if an itemset stays in the top-k largest tiles for more than w consecutive windows and its area is at least twice the area of the k-th tile, then the aTile algorithm finds it. An important point of this theorem is that it does not depends on the value of L. Therefore, even if the probability of false negative or false positive is higher with a small L, the algorithm is still able to mine these tiles that we call the long lasting tiles.

Theorem 4 Let z be a time and J be an itemset. If the area of J in every windows W t with z ≤ t < z +w is larger than two time the area of the k-th largest tile in this window (i.e., A t (J) ≥ 2A t (I t k )), then there is a time z-w < t * < z+w such that J (or one of its supersets) is in C t * and is never pruned by aTile.

Proof. We define hypothesis H as: t * does not exist. We will show that if H is true, it leads to a contradiction and thus the theorem must be true. We denote by [x, y] the subsequence of transactions T x • • • T y . The occurrence times of J in the subsequence [z-w+1, z+w-1] are denoted by z-w < t 1 < . . . < t k < z+w.

If H is true, then for every t i , itemset J and its supersets must be pruned from C ti at some time t i no later than t i + w -1, i.e., t i < t i + w.

Let S i be the subsequence

[t i , t i ] = T ti • • • T t i . If for all z -w < t i ≤ z, we have t i ≤ z then we define W 0 = W z . Otherwise, we take W 0 = W t max where t max = max{t i | t i ≤ z}.
We now construct a set S of these subsequences S i such that every subsequence of S is included in W 0 , every occurrence of J is in at least one of these subsequences, and at most two of these subsequences intersect at any given time. We start with S = ∅ if W 0 = W z or with S = {[t max , t max ]} otherwise. We scan the window W 0 from left to right. If there is an occurrence of J at time t i not already in a subsequence of S, then we add S i in S. All the added subsequences are disjoint by construction, only the last added one may intersect with (t max , t max ).

For every S i ∈ S, let A i be the area of the k-th largest tile of S i used to prune J, i.e., J * < A i . Since by Eq. 2 J * ≥ wA Si (J)/(t i -t i + 1), we have wA Si (J)/(t i -t i + 1) < A i . If we define A m = max A i , then A Si (J) < A m (t it i + 1)/w and by summation for all S i ∈ S: i A Si (J) < A m i (t i -t i + 1)/w. Since all occurrences of J are covered by at least one S i , i A Si (J) > A W0 (J) and since at most two subsequences S i from S intersect at any given time, i (t i -t i + 1) ≤ 2w and thus A W0 (J) < 2A m . Since A m is the size of the k-th tile of one of the S i ⊆ W 0 , it is less than the size of the k-th tile in the whole window W 0 . Finally, A W0 (J) is strictly less than two time the area of the k-th tile in W 0 which is a contradiction.

Experiments

In this section, we perform experiments with two real-life data streams to compare the proposed algorithms to the baseline approaches with respect to the efficiency and the accuracy of the results. The two datasets are: The datasets and the source codes of aTile and cTile in C++ are available for download 5 . The baseline algorithms for comparison are as follows:

-Tile: the original implementation of the tiling algorithm for static database [START_REF] Geerts | Tiling databases[END_REF]. In order to adopt this algorithm for a data stream, we use T ile to recalculate the top-k largest tiles from scratch whenever the window is sliding. -sTile: a sampling based technique proposed in [START_REF] Boley | Direct local pattern sampling by efficient two-step random procedures[END_REF]. The sT ile algorithm samples N itemsets from the window such that each itemset is sampled with probability proportional to the area of the itemset. The top-k largest tiles are extracted from the samples every time the window is sliding.

Accuracy

We performed experiments to demonstrate the effectiveness of the aTile algorithm in term of accuracy. In the first experiment the window size is set to w = 5000 for the Kosarak dataset and w = 10000 for the AOL dataset. Figure 4 shows the average accuracy calculated as the precision of the top-k largest tiles returned by the aTile algorithm when L is varied. Even for very small L, the accuracy is very high, e.g. 99% in the Kosarak dataset and 93% in the AOL dataset. The accuracy increases and reaches 100% accuracy when L is increased.

In the same figure, we can also see that the results are similar when k is varied.

In order to compare to the baseline algorithm sTile, we plot the accuracy (y-axis) versus average update time (x-axis) in Figure 3. We varied both L for the aTile algorithm and the number of samples of the sTile algorithm. Recall that the sTile algorithm depends on the number of samples it collects from the window. When the number of samples increases, so do the accuracy and the update time. This fact is illustrated in Figure 3 in which the sTile algorithm is significantly less accurate than the aTile algorithm given about the same average update time.

The accuracy of sTile may be negatively influenced by the fact that sTile has a very low probability to find each top-k tile. When we calculated the average sum of the area of the top-10 tiles in the Kosarak dataset for w = 5000, we observed that it was significantly lower for sTile than for the aTile algorithm (Figure 5). This confirmed that sTile is not able to find the large tiles with a reasonable number of samples.

Efficiency

Figure 6 shows average update time of four algorithms for k = 10 and different values of w. The number of samples in the sTile algorithm was set to N = 80000 (when N is larger, the algorithm become significantly slower). We set L = w 10 (as shown in Sect. 6.1, the accuracy is then always above 85%). In term of update time, aTile is up to an order of magnitude (50x) faster than the cTile algorithm and about two to three orders of magnitude faster than the sTile (100x) and the Tile (1000x) algorithms. The speed-up increases with the window size.

In Figure 7, the last plot shows the number of candidates kept in the summary of the aTile and the cTile algorithm for varying values of w (the result of the AOL dataset is omitted because it is very similar to the results of the Kosarak dataset in this experiment). The aTile algorithm is more memory efficient than the cTile algorithm as the number of candidates it keeps in the summary is much lower.

Finally, in Figure 7 we show the running time when the aTile, cTile and Tile algorithms are used to find the top-10 largest tiles in a static corpus with varying size. The purpose of the experiment is to see whether the aTile and the cTile algorithm can find the large tiles in a static corpus more efficiently than the Tile algorithm. For the AOL dataset, when most of transactions are small, the Tile algorithm is the fastest one. However, for the Kosarak dataset, when the average transaction size is larger, aTile outperforms both Tile and cTile. Therefore, the aTile algorithm can not only be used for mining large tiles from a data stream but also to mine large tiles from a database efficiently.

Application to Topics Monitoring in the AOL Query Stream

In order to show a potential application of the work we created a demo video6 to visualize the top-k largest tiles of the AOL query stream. Each snapshot of the video corresponds to a list of the largest tiles extracted from a sliding window. For example, in the first snapshot the keyword "real estate" is the most important tile. However, it becomes less important in the second and the third snapshot. Meanwhile, new important tiles such as "high school" is emerging as search for the "high school musical movies" increases in 2006. With this demo, users can track the dynamic of important search topics online.

Application to Tracking in Videos

We investigate how mining top-k largest tiles in a data stream can be useful for analysis of videos and, in particular, for tracking. A more detailed description of this application can be found in [START_REF] Calders | Analysis of videos using tile mining[END_REF]. We worked on a real video made of 5619 frames. This video is shot from a car while following another car (the main object). The main object is present in almost all the frames of the video. If we are able to represent the video as a data stream and use our tile mining algorithm on it, it should discover the car as a large tile or a set of large tiles.

We used the segmentation algorithm7 [START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF] to generate a stream of graphs (one graph per frame). Each graph is a Region Adjacency Graph (RAG) where each node is a region (a set of adjacent pixels) and two nodes are connected if their regions are adjacent. This algorithm performs a temporal segmentation, which means that a given region (and thus, the corresponding node) is present in several successive frames. However, a single region is not enough to track an object in the video: due to change in the object pose or illumination, some regions will split, or merge or disappear and so do the corresponding nodes. We then build a transaction consisting in the set of nodes of each RAG.

Mining tiles on this data stream results in tiles containing regions spread all over the frames. Moreover, the top-k tiles are generally very similar (only differing by one or two regions). Indeed, by removing the information about the adjacency of the regions, it was not possible to find meaningful tiles. The main problem is that the regions composing the car (or more generally, an object to be tracked) must be "close" to each other and this was not taken into account in the original setting. Therefore, we added a spatial constraint to the problem. We chose to mine tiles that are star subgraphs of the RAGs. A star subgraph is a graph with nodes i 0 , i 1 , . . . i n where i 1 , . . . , i n are nodes adjacent to i 0 in the RAG. This formalization ensures that the items in a tile are adjacent to a "center" item i 0 . Moreover, it is easy to integrate this constraint into our algorithms: the intersection of two star graphs with the same center is just the intersection of their set of nodes. The algorithm is otherwise unchanged. We extracted the top-30 tiles with a window of 300 frames (see an example on Fig 10). Then we checked if the extracted tiles could be useful to track the main object (car). We manually drawn the bounding box containing the car in each fifth frame. Then, we selected, for each frame, the tile with the best precision. We end up with a set of 69 tiles with an overall precision of 0.92 and recall of 0.74 on the whole video (and even a recall of 1 on the 2000 first frames). This means that the tiles can actually be used as a high level feature for tracking. Of course, in a real application, the bounding box is not known. But techniques like those described in [START_REF] Diot | Graph mining for object tracking in videos[END_REF] can be used to build tracks from the extracted tiles.

Conclusions and Future Works

In this paper, we proposed two algorithms for mining the top-k largest tiles from a data stream with a sliding window model. The first candidate-based algorithm cTile solves the problem exactly but the update time becomes important when the window size increases. The second one is an approximation algorithm with theoretical bounds on the error rate. Experiments with two real-life datasets show that the approximation algorithm can find large tiles with high accuracy while being an order of magnitude faster than the candidate based algorithm and two to three order of magnitude faster than the other baseline algorithms.

We also show potential applications of large tiles mining in monitoring emerging popular topics in a search engine query log or in a tracking problem. A possible extension for future work is to consider mining different types of tiles with constraints for meaningful applications.

Fig. 2 .

 2 Fig. 2. An example of the summary maintained by the cTile algorithm for sliding windows with size w = 5. The summary contains candidate itemsets which can become a top-k tile in the future. Every candidate is associated with a lower-bound and an upper-bound on the area which are used to prune the candidate list.

Figure 2 (

 2 upper-part) shows an example of a sliding window with size w = 5. Each candidate J in a candidate list C i is associated with two numbers corresponding to the lower-bound and the upper-bound on the area of the candidate. For instance, for the candidate J = ac associated with the candidate list C 3 the lower-bound is J -= 4 because the area of ac in the set of currently observed transactions T 3 , • • • , T 5 is 2 * 2 = 4. For any future window in which the transaction T 3 is not expired, the area of ac is at least as large as 4. Meanwhile J + = 8 because the area of ac in any future window that contains T 3 is at most 4 + 2 * 2 = 8. The area is equal to J + only when the new transactions T 6 and T 7 both contain ac.

Figure 2 (

 2 bottom part) shows one update step of Algorithm 1 when a new transaction T 6 = {a, c} is appended to the window and T 1 together with C 1 are removed from the summary. First, T 6 is intersected with all the existing candidates in C 2 , C 3 , • • • , C 5 to add new closed sets to the candidate lists.

Algorithm 1 1 : 2 :

 112 cTile(S t ) Input: A stream of transactions T1, T2, • • • , T∞, a sliding window size w and a parameter k Output: Summary C for calculating top-k largest tiles 3: C ←-{C1, C2, • • • , Cw} //candidate lists 4: for t = w → ∞ do 5:

Fig. 3 .

 3 Fig. 3. Accuracy versus average update time in seconds of the aTile algorithm and the sTile algorithm when L and the number samples increases

Fig. 4 .

 4 Fig.4. The average accuracy of the aTile algorithm when L is varied. In the table: the average accuracy of the top-10 largest tiles returned by aTile when the window size is varied and L = 0.1w

Fig. 5 .Fig. 6 .

 56 Fig. 5. The average sum (larger is better) of the area of the top-10 tiles returned by the aTile and the sTile algorithms.

Fig. 7 .

 7 Fig.7. Running time according to the database size when the algorithm is used to mine large tiles in a database instead of a stream. The right most subplot shows the number of candidates in the summary when the window size increases.

Fig. 8 .

 8 Fig.8. Top-10 largest tiles in windows (ordered by timestamps) from the AOL query stream. Larger words correspond to tiles with bigger area. The tile "high school" emerges in the second window related to searches for the "high school musical movies".

Figure 8

 8 Figure8shows example snapshots taken from the demo video. Three snapshots of the demo (ordered by timestamps) show the evolution of the largest tiles overtime. Each snapshot is visualized by the wordcloud tool in R. Larger words correspond to tiles with larger area.For example, in the first snapshot the keyword "real estate" is the most important tile. However, it becomes less important in the second and the third snapshot. Meanwhile, new important tiles such as "high school" is emerging as search for the "high school musical movies" increases in 2006. With this demo, users can track the dynamic of important search topics online.

Fig. 9 .

 9 Fig. 9. Frame of the video. Fig. 10. top-30 star-shaped tiles in the segmented frame.

www.twitter.com

http://www.win.tue.nl/ ~lamthuy/tile.htm

http://www.youtube.com/watch?v=3UCjs9d91_g

http://www.videosegmentation.com

Acknowledgements

This work is supported by the Netherlands Organization for Scientific Research (NWO) through the project Mining Complex Patterns in Stream (COMPASS) and the French ANR project SoLStiCe (ANR-13-BS02-0002-01). We would like to thank Dr. Bart Goethals for providing us the implementation of the Tile algorithm and Pr. Toon Calders for his precious comments.


Tri Van

Didier Nguyen

Gildas Georges

Besancon

Van Tri Nguyen
email: van-tri.nguyen@gipsa-lab.grenoble-inp.fr

Didier Georges
email: didier.georges@gipsa-lab.grenoble-inp.fr

Gildas Besancon
email: gildas.besancon@gipsa-lab.grenoble-inp.fr

State and parameter estimation in 1-D hyperbolic PDEs based on an adjoint method

Keywords: State estimation, Parameter estimation, Hyperbolic system, Adjoint method, Inverse problem, Saint-Venant model, Lighthill-Whitham-Richards model

Introduction

In nature and ordinary life, one can find a lot of physical laws described by hyperbolic partial differential equations of order one such as water flow, traffic flow, gas dynamics or electrical lines for instance. Researchers in control have investigated a large number of problems concerning this type of systems with different purposes for example to design an infinite-dimensional nonlinear predictive control for open-channel flow [START_REF] Georges | Infinite-dimensional nonlinear predictive control design for open-channel hydraulic systems[END_REF]; to model and control a dam-river system [START_REF] Litrico | Robust continuous-time and discrete-time flow control of a dam-river system. (I) Modelling[END_REF], [START_REF] Litrico | Robust continuoustime and discrete-time flow control of a dam-river system. (II) Controller design[END_REF]; to investigate the receding horizon boundary control applied on Lighthill-Whitham-Richards traffic flow model to avoid shock waves [START_REF] Thang | Receding horizon boundary control of nonlinear conservation laws with shock avoidance[END_REF]; to optimally control traffic in highway network using linear programming [START_REF] Li | Optimal traffic control in highway transportation networks using linear programming[END_REF]; to apply the predictive control method on gas jet flames distribution [START_REF] Sun | Dynamic System Control of Gas Jet Flames Distribution by Predictive Control Method[END_REF]. One of the important problems arising in the simulation and control of hyperbolic systems is that uncertainty of the initial condition and empirical parameters can cause large errors and inconsistency between the output of control system and the real one. This motivated few studies on observer design, and even output feedback control, going back to works of Christofides et al [START_REF] Panagiotis | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Panagiotis | Distributed output feedback control of two-time-scale hyperbolic PDE systems[END_REF] for instance, with pole placement and Kalman designs , up to [START_REF] Hasan | Adaptive Boundary Control and Observer of Linear Hyperbolic Systems With Application to Managed Pressure Drilling[END_REF] more recently, with backstepping approaches (as also in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2x2 linear hyperbolic system[END_REF]). In the present paper, the purpose is to develop a method based on optimal control theory to optimally estimate the initial condition and distributed parameters in such 1-D hyperbolic systems. Some authors studied constant parameter estimation in hyperbolic system such as Becker et al. who used the method of influence coefficient to estimate Manning roughness coefficient in an unsteady open-channel flow [START_REF] Becker | Identification of parameters in unsteady open channel flows[END_REF]; H. Longxi investigated a complex method to estimate the values of all roughness of a channel network [START_REF] Longxi | Parameter estimation in channel network flow simulation[END_REF]; Y. Ding et al. proposed an adjoint analysis method to find out roughness coefficient in shallow water [START_REF] Ding | Identification of Manning's Roughness Coefficients in Shallow Water Flows[END_REF]. In the case of distributed parameters, Y. Ding et al. also considered the same approach and applied it on the multi-reaches chan-nel flow network [START_REF] Ding | Identification of Manning's roughness coefficients in channel network using adjoint analysis[END_REF]; Richard et al. considered a numerical scheme used to solve parameter identification issue for 1-D hyperbolic system [START_REF] Ewing | A direct method for parameter estimation in a hyperbolic partial differential equation[END_REF]; whereas Wenhuan investigated a quasi-Newton method to deal with the same problem [START_REF] Yu | A Quasi-Newton Method for Estimating the Parameter in a Nonlinear Hyperbolic System[END_REF]. In a former paper of Bagchi et al. [START_REF] Bagchi | Simultaneous ML estimation of state and parameters for hyperbolic systems with noisy boundary conditions[END_REF], the parameter and system state of a discrete-time hyperbolic system with noisy boundary condition are estimated simultaneously based on the maximization of a likelihood function. More recently, we proposed an adjoint method to estimate the initial conditions in an overland flow described by a one dimensional Saint-Venant equation [START_REF] Van Tri Nguyen | Optimal state estimation in an overland flow model using the adjoint method[END_REF]. Notice that contemporaneously with it, Hasan et al. investigated a moving horizon technique to estimate state and constant parameter in a 2x2 linear hyperbolic system based on a distributed model for drilling application [START_REF] Hasan | Moving horizon estimation in managed pressure drilling using distributed models[END_REF]. In that study, the model was discretized first and the adjoint method was applied to the resulting finite dimensional system. In the present paper, which can be considered as an extended version of [START_REF] Van Tri Nguyen | Optimal state estimation in an overland flow model using the adjoint method[END_REF] to general hyperbolic systems, the adjoint method is formulated and solved directly on nonlinear infinite dimensional models. More precisely, we deal with a parameter and state estimation approach in a one-dimensional nonlinear hyperbolic system, with a variable denoted by u(x, t), and a flow denoted by f (u(x, t), x), respectively depending on some initial condition denoted by u i 0 (x), and some distributed parameters denoted by α i (x), both supposed to be unknown, and thus needing to be estimated. The adjoint analysis is formulated with the original infinite dimensional system, to connect the sensitivity of variables needing to be estimated with the system model and the measurements. These points are also the main contributions of the present paper, and to the best of our knowledge, there are very few researches realized with such a spirit. The rest of this paper is organized as follows: section 2 describes the dynamics of system and the formulation of optimal estimation problem. The adjoint method is applied to the optimization problem to get the adjoint system and gradient of estimated variables is presented in section 3. In section 4, two illustrative examples dealing with parameter and state estimation in traffic flow and overland flow system are presented. Some conclusions and perspectives are given at the end of the paper.

Estimation problem statement

System dynamics

Let us consider a general form of 1-D hyperbolic system of variable u(x, t) and flow f (u(x, t), x) a function of u(x, t) and x. The spatial variable x and time variable t belong to the set (x, t) ∈ [0, L] × R + , and the system reads:

       ∂u(x, t) ∂t + ∂f (u(x, t), x) ∂x = g(x, t) u(x, 0) = u i 0 (x) u(0, t) = u b 0 (t) (1)
where the function u b 0 (t) is a predefined boundary condition, the function u i 0 (x) denotes the initial condition, and g(x, t) is a known function. Notice that if it is clear enough in the sequel, notation f and u will be used instead of f (u(x, t), x) and u(x, t) in order to shorten the equation length. Let us assume, without any loss of generality, that the function f (u, x) can be written in the following form for some vectors α = [α 1 (x)...α i (x)...α K (x)] T and ϕ = [ϕ i (u)...ϕ i (u)...ϕ K (u)] T .

f = K i=1 α i (x)ϕ i (u) (2)
On this basis, let us consider the problem of estimating time and space evolution of u when initial condition u i 0 (x) is unknown, together with function parameter α(x). Once the initial state and parameter of system are successfully recovered, all transient state profiles of the system can be fully rebuilt by simulation.

Optimal estimation problem

For the estimation problem of parameter and state in the considered hyperbolic system, one can use two main approaches: empirical procedures or minimization approach based on optimal control theory. The first approach uses direct empirical formula with observation data to get the parameter value and is suitable only for simple parameter estimation. In the present work, due to the complexity and non linearity of estimation problem, we use the second approach to minimize the errors between simulations and some lumped observation values of variable u(x, t). In other words, we minimize a cost function J defined as follows:

J = 1 2 N j=1 T 0 L 0 δ A (x -x j )u dx -u meas j (x j , t) 2 dt + 1 2 ε 1 L 0 u i 0 (x) -u i 0F (x) 2 dx + 1 2 ε 2 K i=1 L 0 α i (x) -α iF (x) 2 dx
(3) where T = optimization horizon (hours); L = considered spacial length where the system takes place (m); N = number of observation values of u(x, t); u meas j (x j , t) = measured value of u(x, t) at observation position x j with x j ∈ [0, L]; α iF (x) = first guessed values of parameters; u i oF = first guessed values of initial condition; ε 1 and ε 2 = weighting factor applied to first guessing term to calibrate the estimated value and the guessed one and adjust the scale of objective function. The term δ A (x -x j ) is an approximation of Delta-Dirac function described by a Gaussian function with a very small variance σ 2 as

δ A (x -x j ) = e -(x-xj ) 2 /σ 2 .
3 Adjoint-based approach

Variational analysis

From the previous analysis, the optimal values of state and parameter must minimize the cost function in equation [START_REF] Panagiotis | Feedback control of hyperbolic PDE systems[END_REF] and satisfy also the system equation ( 1). This constraint and the continuity of the first partial derivative of both system dynamics and cost function lead to use the Lagrange multiplier with Lagrangian variable λ(x, t), which allows to combine the system equation and cost function into only one new cost functional L(λ, ρ, α), shortened to L, as follows.

L = J + T 0 L 0 λ ∂u ∂t + ∂f ∂x -g(x, t) dxdt A (4)
To deal with this optimization problem the common adjoint method is used in order to obtain the adjoint system and establish the gradient of cost functional with respect to the parameters and state needing to be estimated. These gradients describe the sensitivity of cost functions under the constraints of governing equation and variation of parameters. First of all, we take the first variation of the cost functional with respect to the system variable u, initial condition u i 0 (x) and weighting factor α i (x). The first variation δJ is given in equation [START_REF] Ding | Identification of Manning's Roughness Coefficients in Shallow Water Flows[END_REF].

δJ = N j=1 T 0 L 0 δ A (x -x j ) L 0 δ A (x -x j )udx -u meas j δudxdt + ε 2 K i=1 L 0 α i (x) -α iF (x) δα i (x)dx + ε 1 L 0 u i 0 (x) -u i 0F (x) δu i 0 (x)dx (5)
Then using the integration by part technique the first variation of second term, called A, of L can be obtained and presented in equation [START_REF] Ding | Identification of Manning's roughness coefficients in channel network using adjoint analysis[END_REF].

A = L 0 λu T 0 dx - T 0 L 0 u ∂λ ∂t dxdt + T 0 λf L 0 dt - T 0 L 0 f ∂λ ∂x dxdt - T 0 L 0 λg dxdt (6)
The first variation of A, calculated with respect to system variable u(x, t), initial condition u i 0 (x) and param-eters α i (x), is developed as follows:

δA = L 0 λδu T 0 dx - T 0 L 0 ∂λ ∂t δudxdt + T 0 λ ∂f ∂u δu L 0 dt + K i=1 T 0 λ ∂f ∂α i (x) δα i (x) L 0 dt - T 0 L 0 ∂f ∂u ∂λ ∂x δudxdt - K i=1 T 0 L 0 ∂f ∂α i (x) ∂λ ∂x δα i (x)dxdt (7)
All the terms multiplied by δu in equations ( 5) and ( 7) is collected together and set to zero to cancel all the value of δu. This gives the first order of optimality condition or the adjoint system of adjoint variable λ(x, t), described in equation ( 8):

- ∂λ ∂t - ∂f ∂u ∂λ ∂x + N j=1 δ A (x -x j ) × L 0 δ A (x -x j )udx -u meas j = 0 (8)
where

∂f ∂u = K i=1 α i (x) dϕ i (u) du .
The gradient of objective functional with respect to parameter α i (x), formed by selecting all the terms related to δα i (x) in previous equations of first variation, is not identical throughout the considered spacial domain [0, L]. At positions x = 0 and x = L the gradients of α i (x) are presented in equations [START_REF] Hasan | Adaptive Boundary Control and Observer of Linear Hyperbolic Systems With Application to Managed Pressure Drilling[END_REF] and [START_REF] Hasan | Moving horizon estimation in managed pressure drilling using distributed models[END_REF]. At other positions, this gradient is in equation [START_REF] Lax | Systems of conservation laws[END_REF].

∂L ∂α i (0) = - T 0 λ(0, t) ∂f ∂α i (0) dt (9)
∂L ∂α i (L) = T 0 λ(L, t) ∂f ∂α i (L) dt (10)
∂L ∂α i (x) = - T 0 L 0 ∂f ∂α i (x) ∂λ ∂x dxdt + ε 2 L 0 α i (x) -α iF (x) dx (11)
where

∂f ∂α i (x) = ϕ i (u).
By using the same method, the gradient with respect to initial condition u i 0 (x) is obtained and presented in equation [START_REF] Li | Optimal traffic control in highway transportation networks using linear programming[END_REF].

∂L ∂u i 0 (x) = - L 0 λ(x, 0)dx + ε 1 L 0 u i 0 (x) -u i 0F (x) dx (12)
The system dynamics and adjoint equations must be discretized in order to be solved numerically. The considered spacial domain [0, L] is discretized into a set of smaller sections [x j , x j+1 ]. As a result the gradients of initial condition and parameters α i (x) at position x j , ∀x j = j x where 0 < j ≤ M -1 with x = L/M , M being the number of point, become:

∂L ∂u i 0 (x j ) = -λ(x j , 0) + ε 1 u i 0 (x j ) -u i 0F (x j ) (13)
∂L ∂α i (x j ) = - T 0 ∂f ∂α i (x) ∂λ ∂x x=xj dt+ε 2 α i (x j )-α iF (x j ) (14)
According to the optimization condition, the gradient of u i 0 (x) must be zero to minimize the cost function. We can obtain this condition, through optimization process, by setting all the gradient at each discretization point x j to zero. After forming the adjoint system and gradients, all the remaining terms of equations ( 5) and ( 7) are gathered together and set to zero to satisfy the optimization conditions and to get the initial and boundary condition of adjoint variable λ(x, t).

L 0 λ(x, T )δu(x,T )dt + T 0 λ(L, t) ∂f (u(L, t), L) ∂u(L, t) δu(L, t)dt - T 0 λ(0, t) ∂f (u(0, t), 0) ∂u(0, t) δu(0, t)dt = 0 (15)
where

∂f ∂u (L, t) = K i=1 α i (L) ∂ϕ ∂u (L, t).
The boundary condition of system variable u(0, t) is fixed beforehand, as a result its variation δu(0, t) is zero for all values of t. The first and second terms of equation ( 15) also vanish also if we make λ(x, T ) equal to 0 and λ(L, t) must satisfy equation [START_REF] Longxi | Parameter estimation in channel network flow simulation[END_REF].

λ(L, t) K i=1 α i (L) ∂ϕ ∂u (L, t) = 0 ∀(x, t) ∈ [0, L]×R + (16)
They are the necessary conditions to solve the adjoint system [START_REF] Georges | Infinite-dimensional nonlinear predictive control design for open-channel hydraulic systems[END_REF]. We can notice that the boundary condition λ(L, t) of adjoint system depends on the value of the partial derivative of function ϕ(u(L, t)) with respect to system variable at position L. If this derivative does not vanish, then, λ(L, t) must be equal to zero. This makes the gradient of parameter α i (L) in equation [START_REF] Hasan | Moving horizon estimation in managed pressure drilling using distributed models[END_REF] equal to zero and we can not estimated the value of parameter α i (L). We can also notice that the gradient of initial condition u i o (L) in equation ( 13) becomes much smaller compared to other positions because we only have the correction term with calibration factor 1 . With small 1 the estimation of u i o (L) does not converge. However, an ε 1 large enough can reduce the biases on the estimation of u i o (x i ), ∀x j = j x where 0 ≤ j ≤ M -1 with x = L/M , M being the number of point. In two numerical examples presented later on in section 4, the problem of null value of λ(L, t) is met. A temporary solution is to consider that the values of u i o (L) and α i (L) are measured, and to only deal with the estimation distributed in the internal spatial section ∀x j = j x where 0 ≤ j ≤ M -1 with x = L/M , M being the number of point. A specific analysis will be carried on to investigate more deeply this incompleteness of adjoint method for estimating the value of variable at the end of the considered domain.

Numerical implementation

Explicit Lax-Wendroff scheme for discretization

The Lax-Wendroff scheme belongs to the class of conservative schemes and was presented in 1960 by P.D. Lax and B. Wendroff [START_REF] Lax | Systems of conservation laws[END_REF]. It is a numerical method for solving the systems of one dimension hyperbolic conservation laws with second order of accuracy in both space and time. It can be easily extended to the case of general 1-D hyperbolic system (1) with non homogeneous term g(x, t). The formula of this scheme can be derived by a variety of approaches including Taylor series expansion in time of variable u(x, t):

u(x, t+ t) = u(x, t)+ t ∂u ∂t + t 2 2 ∂ 2 u ∂t 2 +O( t 2 ) (17)
The first and second partial derivative of u(x, t) with respect to t can be replaced by space derivatives of and others term of equation ( 1), which are assumed to exist with enough smoothness. The value of variable u at time t + t is as follows, neglecting the higher order terms O( t 2 ) in Taylor series expansion :

u(x, t + t) = u(x, t) + t - ∂f ∂x + g + t 2 2 × - ∂ ∂f ∂u ∂x - ∂f ∂x + g - ∂f ∂u - ∂ 2 f ∂x 2 + ∂g ∂x + ∂g ∂t (18)
All the space and time derivatives in equation ( 18) are approximated by central differences. By denoting u(x + i x, t+n t) by u(i, j) and ∂f ∂u = h, the Lax-Wendroff scheme for discretizing equation (1) becomes:

u(i, n + 1) = u(i, n) + t - f (i + 1, n) -f (i -1, n) 2 x + g(i, n) + t 2 2 - h(i + 1, n) -h(i -1, n) 2 x × - f (i + 1, n) -f (i -1, n) 2 x + g(i, n) -h(i, n) × - f (i + 1, n) -2f (i, n) + f (i -1, n) x 2 (19)
+ g(i + 1, n) -g(i -1, n) 2 x + g(i, n + 1) -g(i, n -1) 2 t (20)
The conventional truncation error analysis shows that this 3-point scheme is second-order accurate in space and time. Moreover, its stability and convergence is guaranteed by the Courant-Friedrichs-Lewy (CFL) condition.

In our case of one-dimensional hyperbolic system, this

CFL reads | u max | t x C max .
The typical value of coefficient C max for this explicit Lax-Wendroff scheme is 1. It is worth noting that this scheme has also been chosen after comparison of its computational burdern with other possible approaches. In particular, an alternative implicit four-point Preissman scheme was considered (see [START_REF] Van Tri Nguyen | Optimal state estimation in an overland flow model using the adjoint method[END_REF]) but appeared to be more time consuming. For example, to numerically simulate the traffic flow of subsequent subsection 4.1.1 (with 10 space discretization steps and 300 time discretization steps), it takes 0.24 seconds with Lax-Wendroff scheme, and 10.67 seconds with Preissmann scheme, using the same simulation parameters and same computation power (Intel(R) Core(TM)2 Duo CPU P8600 @2.40 GHz, 3Go RAM with Matlab R2011b 32 bit).

Method for optimal estimation problem resolution

The optimization problem is solved by using optimization tool called f mincon in Matlab. The gradient of cost functional with respect to the initial condition and parameters got from equations ( 9), ( 10), ( 14) and ( 13) is provided to this function as well. In order to prevent the estimated variables for u i 0 (x) and α i (x) to become negative and be out of the limits describing the physical meaning, we imposed also a linear constraint. It is worth to notice that the f mincon tool has itself a finite difference method to calculate the gradient of estimated variable, which means that it is not obligatory to provide it via adjoint method in some simple estimation problem. In our two examples, due to the estimation size and the complexity of considered nonlinear hyperbolic models, the optimization problem converges with long estimation time (in the case of traffic flow) or does not converge (in the case of overland flow) without feeding gradients.

Illustrative examples

In order to validate the effectiveness of the proposed method, we present hereafter two numerical examples taken from two different application fields: traffic flow dynamics on the one hand and overland flow dynamics on the other hand. The validation procedure is developed using free-error approach which means that all the measurements or observed values of system variable u(x, t) is got from the exact simulation of system with optimal initial state and parameters without noise adding. These measurements are used as input of optimization process.

The accuracy of estimation result depends on the number and position of observation values. In this work, the number of sensors and their positions are selected manually, and a full sensitivity analysis of distributed state and parameter estimation with respect to the observations is left for future works. The traffic flow here-considered is described by the well-known Lighthill-Whitham-Richards (LWR) model [START_REF] Michael | On kinematic waves. ii. a theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF] which was firstly introduced by Lighthill and Whitham in 1955 and complemented by Richards in 1956. This is a macroscopic traffic flow model in which the movement of an individual vehicle is neglected, and only a set of vehicles with its characteristics like density and average velocity, is considered to formulate the dynamics of traffic flow, as presented in equation [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2x2 linear hyperbolic system[END_REF].

∂ρ(x, t) ∂t + ∂ ρ(x, t)v(x, t) ∂x = 0 (21)
where x ∈ R; t ≥ 0 are space and time variables; ρ(x, t) is the density of vehicles (vehicle/km) in x at time t; v(x, t) is the average velocity of vehicles in x at time t (km/hour). We make an assumption that the considered highway is long enough to be divided into a set of small sections, each of them with its own parameters. Moreover, on each section the velocity only depends on density ρ(x, t) through the simple function [START_REF] Sun | Dynamic System Control of Gas Jet Flames Distribution by Predictive Control Method[END_REF] below. This means that when the highway is empty or ρ = 0, the vehicles can get the maximal velocity v = v max . When the highway is full of vehicles ρ = ρ max , the vehicles can not move v = 0.

v(x, t) = v max 1 - ρ(x, t) ρ max (22)
The parameters v max and ρ max are supposed to be constant. This assumption is reasonable in real high way whose characteristics do not vary a lot along its length. As a result, the flow function is f (ρ, x) = v max ρ(x, t)v max ρ 2 (x, t)/ρ max . We consequently have the vectors

of parameters α = [v max v max /ρ max ] T and ϕ(ρ) = [ρ(x, t) -ρ 2 (x, t)] T and for f (ρ(x, t), x): ∂f (ρ(x, t), x) ∂α 1 = ρ(x, t) ∂f (ρ(x, t), x) ∂α 2 = -ρ 2 (x, t) ∂f (ρ(x, t), x) ∂ρ(x, t) = α 1 -2α 2 ρ(x, t) (23)
Let us insert the partial derivative of flow function ∂f (ρ(x, t), x) ∂ρ = α 1 -2α 2 ρ(x, t) into equation ( 8). Then we can get the adjoint system corresponding to the direct traffic flow in equation [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2x2 linear hyperbolic system[END_REF]. The gradient of cost functional L with respect to initial condition of vehicle density is the same as equation [START_REF] Michael | On kinematic waves. ii. a theory of traffic flow on long crowded roads[END_REF]. The sensitivity of constant parameters α 1 and α 2 is a little bit different with the previous general calculation and described as follows.

∂L ∂α i = T 0 λ ∂f ∂α i L 0 dt - T 0 L 0 ∂f ∂α i ∂λ ∂x dxdt + ε 2 L 0 α i -α iF dx (24)
with the derivatives ∂f ∂α 1 = ρ(x, t) and ∂f ∂α 2 = -ρ 2 (x, t)

Estimation result

The estimation process for estimating the optimal values of initial condition of vehicles density and two parameters α 1 and α 2 is realized on a straight traffic road with characteristics presented in Table A

(x) = 20 (vehicle/km) ∀ x ∈ [0...L], α 1 = v max (x) = 100 (km/hour), α 2 = v max /ρ max = 0.25 (km/hour × km/vehicle).
Both the first guessed value of initial condition and two parameters is set to zero. The stopping criterion of optimization tool f mincon is whenever the tolerance on the cost function value J or tolerance on estimation variable is smaller than the termination tolerance values. We took just 1 observation value at the end of considered domain, from a simulation of direct system [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2x2 linear hyperbolic system[END_REF] with real value of initial state and parameters. One can observe in Fig. 1 the estimation results for ρ(x, 0). It converges to the respective real values indeed with some small bias values at some positions in the considered spatial section ∀x j = j x where 0 ≤ j ≤ M -1 with x = L/M , M being the number of point. The estimated value of parameters α 1 and α 2 are 149.8782 and 0.4988 respectively. These estimation biases are possibly caused by the numerical errors on direct and adjoint system simulation and convergence limitation of BFGS method. The estimated value of maximal vehicles density can be easily found by simple division ρ max = α 1 /α 2 . At the end of optimal estimation process the found minimal value of cost function J is 0.0074 slightly larger than zero, the theoretical value, due to estimation biases. The relative errors between the real values of estimation variables and estimated ones are 0.2330%, 0.3011%, and 0.0812% respectively for v max , v max /ρ max and ρ(x, 0) which are pretty small. After some simulation realized to determine the sensitivity of estimation result with respect to the number and position of measurements, x(km) we can conclude that the more measurements we have, the more accurate result we obtain. But the accuracy improvement is not noticeable while being compared to the cost we pay for sensor number increment. Moreover, in this simulation, if the only one sensor is placed on the last position of considered domain, we obtain better optimization result than other positions. By definition, an overland flow is water, caused by excess rainfall or surface runoff, that flows over ground surface and towards some hydraulic channels. The dynamics of this type of water flow can be derived from the well-known Saint-Venant equations by some assumptions about small bottom slope, local and convective acceleration [START_REF] Van Tri Nguyen | Optimal state estimation in an overland flow model using the adjoint method[END_REF]. Finally the dynamics of overland flow is as follow:

ρ i o (x)(vehicle/km) Estimated ρ i o (x) Real ρ i o (x) ρ i o (x) initial
       ∂h ∂t + ∂f (h, x) ∂x = r(x, t) -i(x, t) h(x, 0) = h i 0 (x) h(0, t) = h b 0 (t) (25)
where x = distance, (m); t = time, (s); h = water flow depth, (m); f (h, x) = h 5/3 S 1/2 0 /n(x) = the flow per unit width (m 2 /s) with the spacial distributed Manning roughness coefficient n(x); r = rainfall intensity, (m/s); i = infiltration intensity, (m/s); S 0 = bed slope, (m/m). The infiltration i(x, t) [START_REF] Mays | Water resources engineering[END_REF], caused by variable rainfall, can be described by the famous simple Green-Ampt model in equation ( 26

):    i(t) = K i Ψη(1 -θ) I(t) + 1 if t > t p i(t) = r if t ≤ t p (26)
where K i = effective hydraulic conductivity of the soil of flow, ( cm s ); Ψ = soil suction at wetting front, (cm); η = soil porosity, (%); θ = relative initial soil moisture, (unitless); I(t) = the cumulative infiltration, the accumulated depth of water infiltrating during the time period t, (cm); t p = ponding time, (s) calculated by

t p = ΨK i η(1 -θ) r(r -K i )
. By comparing this to the general form of flow in equation ( 2), we have the vector of only one parameter α(x) = S 1/2 0 /n(x) and the corresponding function of system variable ϕ(h) = h 5/3 . By inserting the derivative ∂f (u, x) ∂h = 5 3 α(x)h 2/3 into equation [START_REF] Georges | Infinite-dimensional nonlinear predictive control design for open-channel hydraulic systems[END_REF] we can obtain the adjoint system of overland flow.

Replacing the derivative ∂f (h, x) ∂α(x) = h 5/3 into equations ( 9), ( 10) and ( 11), we can get the gradient of cost functional L with respect to this parameter.

Estimation result

This example is assumed to be taking place on some hills with characteristics described in Table . A.2. We need to estimate the initial condition of flow height h(x, 0) ∀x ∈ [0...L] whose real value is depicted in Fig. 2 and distributed parameter α whose real value α(x) ∀x ∈ [0...L], found via equation α(x) = S 1/2 0 /n(x) with chosen values of n(x), is in Fig. 3 . With the same estimation method as before, 2 observation values are taken out of 10 discretized sections (1 at the begin and 1 at the end) of considered domain and used as inputs for optimization process. The initial values of estimated state and parameter are 0.35(m) ∀x ∈ [0...L] and 7 ∀x ∈ [0...L] for, respectively, h(x, 0) and α(x). After 192 iterations, the estimation of h(x, 0) indeed converges to the desired value in Fig. 2 and parameter α(x) also converges in Fig. 3. The final value of cost function J is 1.1478 × 10 -18 which is very small. Taking a look at the average relative errors of estimated values of h i o (x) and α(x), respectively equal to 3.4240 × 10 -8 % and 2.4165 × 10 -8 %, we can see that the estimation accuracy is pretty good. The estimated Manning coefficient n(x), obtained via optimized parameter α(x), has also a small relative error (3.4240 × 10 -8 %).

Conclusions

A general approach for parameter and state estimation in 1-D hyperbolic system based on optimal control theory and adjoint method has been developed and validated in this paper. First of all, the system dynamics and optimal estimation problem from lumped observation values have been formulated. Then the adjoint based method has been introduced to get the adjoint equation and gradients of objective functional, obtained by Lagrange multiplier method. From these gradients, the BFGS method has been used to find the optimal value of all parameter and state of variables which minimize the objective functional. The partial differential equations of both direct and adjoint system have been discretized and solved numerically by means of implicit Lax-Wendroff scheme. Finally, two numerical examples of hyperbolic system, overland flow modeled by Saint-Venant equation and traffic flow modeled by LWR equation, have been provided to illustrate the effectiveness of the proposed approach. Because of its generality, this method can be applied to the state and parameter identification in any 1-D hyperbolic system. The observation data in this work are obtained from the simulation of system with optimal parameters and it my not represent the complexity and uncertainty of natural problem. As a result, a validation with real data or laboratory experiments is in progress to evaluate the applicability of the presented method to real problem. Moreover, more gen-

4. 1

 1 State and lumped parameter estimation for traffic flow 4.1.1 Dynamics of traffic flow system

Fig. 1 .

 1 Fig. 1. Estimation of initial condition ρ(x, 0)

4. 2 4 . 2 . 1

 2421 State and distributed parameter estimation for hydraulic flow Dynamics of over land flow system

Fig. 2 .Fig. 3 .

 23 Fig. 2. Estimation of hi 0 (x)

  .1. The real value of maximal vehicles velocity α 1 (x) is 150 (km/hour) and the parameter α 2 is 0.5 (km/hour × km/vehicle). The high way length L = 100 (km) is divided into 10 sections with spatial discretization step x = 10 (km). The real value of initial condition and parameters can be also found in table A.1. The boundary condition is a sum of four sinusoidal signals with different frequencies and phases. The initial value feed to the optimization tool is ρ init

eral problems such as output-feedback stabilization or optimal control with output constraints and time varying parameters will be part of our future developments.

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) and the MEPIERA project, Grenoble Institute of Technology

A Tables of simulations parameters

Parameter


Haoran Zhang
email: haoranzhang0715@163.com

Xuyang Wei
email: bweixuyang321@163.com

Tengfei Zou

Zhongliang Li

Guocai Yang

Agriculture Big Data: Research status , challenges and countermeasures

Keywords: agricultural big data, cloud computing, challenge, countermeasure

Agriculture data type and amount in human society is growing in an amazing speed which is caused by the emerging of agricultural Internet of things. The growth of agricultural data volume brought many difficulties in storage and analysis. The realization of the big data and cloud computing technologies provides the solution of these problems. The cloud computing and big data technologies can be applied in agriculture to solve the problems in storage and analysis. This paper elaborates the related concept of agricultural big data briefly, and emphasizes some current researches, the challenges that agricultural big data should face in the future and its countermeasures.

Introduction

Until now, the agricultural big data has not been clearly defined. We can reach the definition with two parts -big data and agricultural informationization. The accurate definition of big data also has some disputes. The well-known consulting company McKinsey argues that big data cannot be obtained, stored and managed in a certain time using a traditional database software tools [START_REF] James | Big data: The next frontier for innovation, competition, and productivity[END_REF]. The definition given by wikipedia shows big data is a blanket term for any collection of data sets which are so large and complex and it becomes difficult to process using onhand database management tools or traditional data processing applications [2]. there is another view that big data should meet three terms which are Variety, Velocity and Volume [START_REF] Grobelnik | Big Data Tutorial[END_REF]. Some people has a different view. They think that three terms can't describe the characteristics of big data, so they come up with four terms. But there is no unified statement on the fouth characteristic. The International Data Corporation IDC believes big data should also has a term of Value [START_REF] Barwick | The" four Vs" of Big Data[END_REF]. The IBM corporation argues that big data should be considered with Veracity [START_REF]What is big data[END_REF]. I think that the term of value makes the definition of big data more accurate. Although the value of big data is sparse, there is great potential inherent value. Agricultural informationization is a dynamic concept. It improves the comprehensive productivity of agriculture and the overall productivity of agriculture using modern information technology and information systems [START_REF] Zeng | Effective Methods for Agricultural Informationization[END_REF]. In brief, agricultural informationization makes full use of information technology methods, means and process to achieve the goals. Agricultural big data uses concepts, techniques and methods of big data and cloud computing to handle a mount of agricultural data. We can obtain the information which is useful to guide agriculture. Wen Fujiang shows that agricultural big data is related to agricultural production in all aspects. It is a multi-professional data mining process [START_REF] Wen | Strategic significance of data and collaboration mechanisms agricultural research[END_REF]. The implementation of agricultural big data technology is a very important component of agricultural informationization. And it also provides new methods and ideas for agricultural research and agricultural business development. Professor Wang refers that we have many datasets of agriculture. As a result, there is a good deal of agricultural data. But the data standards are not uniform and standardized. How to give a reliable and professional decision is an urgent task [START_REF] Wang | Bottleneck of Agricultural Informatization Development in China and the Thinking of Coping Strategies[END_REF]. In order to solve the problem of the development of agricultural informationization, we should establish a national big data center and make efforts to develop cloud computing and to search agricultural big data mining techniques.

Developing Status of Agricultural Big Data

With the rapid rise of the Internet of things and social networks, the amount of data is growing at an unprecedented rate. According to Winter Corporation's survey, the current amount of data is growing at three-time increases every two years [START_REF] Wintercorp | The Large Scale Data Management Experts[END_REF]. Its growth rate exceeds the growth rate of Moore's Law far. The era of big data is coming. "NATURE" magazine published a special issue named "Big Data: Science in the Petabyte Era" [START_REF]Big Data[END_REF]. This paper described the challenges of big data in many aspects,such as Internet technologies, network economics and so on. "SCIENCE" magazine also launched a special issue named "Dealing with Data"[12] in 2011. It illustrated the important scientific value of big data. European Information Society Science and Mathematics Research journal ERCIM News discussed the management and innovative technology issues in their special issue named "Big Data" [START_REF]Big Data[END_REF] in 2012. Gartner Company showed that the development of big data has entered the peak period in 2013 Hype Cycle for Emerging Technologies. We showed it in Fig. 2-1.

Many transnational corporations, such as IBM,ORACLE and FACEBOOK, are making their efforts to develop big data technologies [START_REF] Li | Research status and scientific thinking of big data[J][END_REF]. Academician Professor Li Guojie noted that we should take the redundant data out of raw data in the data explosion era [START_REF] Li | The New Focus of Information Technology Big Data[END_REF]. In December 2013, China Big Data Technology Conference was held and it mainly discussed the technical means and the commercial value of big data applications. Big data has been a intersection of information science, social science, network science, and many other emerging interdisciplinary field. It has became a hot research point. While agricultural big data is still a fresh vocabularies, but the number of data which is produced in the process of agricultural production is far more than small data. Agricultural data generation method has been changing dramatically. With the wisdom agriculture proposed, widely used sensors promote the development of agricultural Internet of things greatly. Agricultural Internet of things generates a lot of agricultural data. We should use big data technologies to analysis the data. The emergence of cloud computing provide a reliable way to deal with these massive data. At this stage, there is few research on agricultural big data. But large number of researchers have realized the value of agricultural big data, and have also been put into analyzing agricultural big data and using agricultural informations. If we can make good use of agricultural big data, it will be not only a great innovation in the history of human agriculture, but also a pioneering work in human history.

Challenges of Agricultural Big Data

In its infancy, big data technology is facing many challenges, such as wide range of heterogeneous data, problem of real-time, data incompleteness, lack of priori knowledge, private issues and so on. Issues which agricultural big data is facing is consistent with the big data technology. But compared with big data, agricultural big data is not so sensitive in security or privacy issues. Agricultural data mining is targeted at using the result to guide agricultural practices. Therefore, agricultural big data is facing the following issues in sum.

Problem of Agricultural Big Data Storage

Agricultural big data has a very different modality. From the source point of view, the data comes from Internet of things of the radio equipment, the agricultural information websites, and a variety of advanced mobile terminals. From the content point of view, it includes not only statistical data, but also basic information on agriculture-related economic entities, investment information, import and export information and GIS coordinate information. Data types also include Structured data, semi-structured data and unstructured data. It will be a problem worth studying to store heterogeneous data and the ability to read and write because of the different treatment of different storage hardware devices. The heterogeneous hardware is also a problem of storing agricultural big data. There will be a very significant performance differences between different machines in the data center. Different hardware devices have different literacy and processing capabilities. Handing equipment will waste a lot of time waitting the slower storage devices. The linear growth of storage devices and servers will not necessarily bring the linear growth of computing power in this case. The "Barrel Effect" restricts the performance of the entire cluster.

Problem of Agricultural Big Data Analysis

Data analysis is the core of the whole process of agricultural big data, because the value of agricultural big data is produced in the big data analysis process. Currently, there are many problems in agriculture such as food security, soil management, pest forecasting and prevention, soil management and agricultural consumption. We can use the analysis of agricultural big data to solve these problems. Raw data comes from the extraction and integration of agricultural information. We can select all or part of these data to do the research. Conventional analytical techniques are not applicable on processing agricultural big data, such as data mining, machine learning, statistical analysis and other techniques. It is mainly discussed in the following parts.

(1)Traditional data mining algorithms, such as machine learning and other areas, are no longer suitable for agricultural big data. On one hand, algorithms which mine a small amount of data can not be directly applied to big data. On the other hand, agricultural big data has its particularity. The accuracy of the algorithm is no longer the main standard. Algorithm needs to strike a balance between timeliness and accuracy of processing in many cases.

(2)The metrics of the quality of the analysis result is also a major challenge. Big data has complex types. It leads to many problems in designing the indicators to algorithm.

Problem of Agricultural Big Data Timeliness

As the time elapsing, the inherent value of data keeps attenuating. Therefore, timeliness must be considered during the analysis of agricultural big data. Untimely data analysis may result in the production of agricultural disasters, Especially in the weather, environmental conditions associated with data analysis. For example, the occurrence of "Low grain price hurts farmers" event is the result of managing the cost of production and other information not in time. Therefore, the characteristic of timeliness is particularly important in agricultural big data.

Countermeasures of Agricultural Big Data

To solve the problem of agricultural big data, we will give the corresponding countermeasures which will give some guidances for the future work.

Storage of Agricultural Big Data

Agricultural big data storage affects not only the efficiency of data analysis, but also the cost of data storage. Therefore, we need to come up with high-efficiency and lowcost data storage. One solution is using distributed storage mode. Distributed file system is the physical storage resources in the file system and is not necessarily directly connected to the local node. But the node is connected via a computer network. Many companies which have to deal with huge amount of data have their own distributed file system, such as Google's GFS(Google File System) [START_REF] Ghemawat | The Google file system [C[END_REF][START_REF] Mckusick | GFS: Evolution on Fast-forward[END_REF], Taobao File System (TFS) [START_REF] Chucai | TFS Introduction[END_REF]. There are some open source distributed file system,such as HDFS [START_REF] Konstantin | The Hadoop Distributed File System[END_REF],NFS [START_REF] Osadzinski | Network File System(NFS)[END_REF][START_REF] Anderson | Serverless Network File Systems[END_REF]. There is another solution. The research can be concerned in searching the technologies about the acquisition and integration of multi-source ans multi-modal data. In addition, high-speed storage and creating index are also important aspects.

The general solution to solve the problem of heterogeneous hardware is using different storage devices in different aspects in heterogeneous hardware environments. The problem will become very complicated when the scale of heterogeneous environment extends to thousands of clusters.

Analysis of Agricultural Big Data

After years of researches and developments, data mining, machine learning, statistical analysis and other information analysis have been proved to have significant effects for small data. These algorithms can be adjusted to accommodate cloud computing system. But it must be noted that we must consider the characteristics of agricultural big data in the adjustment process of these algorithms. Real-time and predictable characteristics must be considered. It will be a hot spot in the coming period of scientific research.

It is important and difficult to evaluate the results of agricultural big data algorithms.

According to the characteristics of agricultural big data, we can use timeliness as a measure standard. Agricultural big data is massive, so we can use the prior knowledge to test the algorithms. It can measure the quality of the algorithms to a certain extent. It can also measure the reliability of the data results.

Timeliness of Agricultural Big Data

Timeliness is a core demand of agricultural big data analysis. A lot of research are also expanded around this demand. There are three methods to ensure timeliness. The first method is using stream processing mode. Although streaming mode is suitable to real-time system, its application field is relatively limited. Streaming application model focuses on real-time statistical system, online monitoring. The second method is batch mode. In recent years, the development of batch model realtime system has became a hot topic and has achieved a lot of achievements. MapReduce programming model which Google company made in 2004 is the most representative batch mode. The third method is using a combination of stream processing and batch mode. The main idea is to use the MapReduce programming model to achieve stream processing.

Fig. 1 .

 1 Fig.1. Emerging Technologies Hype Circle,2013[START_REF] Fenn | Emerging Technologies Hype Cycle for 2013[END_REF] Recently, big data has been applied to different points of healthcare, manufacturing, transportation, and financial sector in China. With the use of agricultural informationization and Internet of things in agriculture, agricultural big data will become another focus of big data research. Agricultural big data is a application of big data in agriculture. Agricultural big data which is related to all aspects of agricultural production is a issue of multi-disciplinary, cross-sectoral analysis and data mining. The first domestic agricultural big data industry technology innovation strategic alliance was established in Shandong Agricultural University in June 2013. Shandong Agricultural University President Wen Fujiang pointed out that in China, the current big data research is just beginning, but agricultural big data research is leading research[START_REF] Wen | Strategic significance of data and collaboration mechanisms agricultural research[END_REF]. Many companies are targeting agricultural big data opportunities. The New York Times reported that in 2012 the Solum company received 17 million dollars fund form Andreessen Horowitz and other companies. Solum company commits to use data analysis techniques to determine the amount of fertilizer in cultivation. It helps farmers improve productivity and reduce costs through the analysis of agricultural big data. Multinational agricultural biotechnology company Monsanto spent 930 million to takeover the insurance accident weather company Climate Corporation[18]. Climate Corporation predicts the weather which may damage agricultural production through analyzing their massive weather data. According to the forecast, farmers can select appropriate agricultural insurance, to reduce the impact of bad weather in agricultural process. The analysis of agricultural big data plays an indispensable role in agriculture.

Corresponding Author Email: paul.g.yang@gmail.com

Acknowledgment

This work was supported by the Key Technologies R & D Program(2012BAD35B08), China.