extfieldsofstudy:
Computer Science
Engineering

text:
Selective Memory Recursive Least Squares: Uniformly Allocated Approximation Capabilities of RBF Neural Networks in Real-Time Learning

When performing real-time learning tasks, the radial basis function neural network (RBFNN) is expected to make full use of the training samples such that its learning accuracy and generalization capability are guaranteed. Since the approximation capability of the RBFNN is finite, training methods with forgetting mechanisms such as the forgetting factor recursive least squares (FFRLS) and stochastic gradient descent (SGD) methods are widely used to maintain the learning ability of the RBFNN to new knowledge. However, with the forgetting mechanisms, some useful knowledge will get lost simply because they are learned a long time ago, which we refer to as the passive knowledge forgetting phenomenon. To address this problem, this paper proposes a real-time training method named selective memory recursive least squares (SMRLS) in which the feature space of the RBFNN is evenly discretized into a finite number of partitions and a synthesized objective function is developed to replace the original objective function of the ordinary recursive least squares (RLS) method. SMRLS is featured with a memorization mechanism that synthesizes the samples within each partition in real-time into representative samples uniformly distributed over the feature space, and thus overcomes the passive knowledge forgetting phenomenon and improves the generalization capability of the learned knowledge. Compared with the SGD or FFRLS methods, SMRLS achieves improved learning performance (learning speed, accuracy and generalization capability), which is demonstrated by corresponding simulation results.

capability), appropriate training methods are essential. For neural network control, the RBFNN is usually employed as an adaptive controller or observer to interact with the plant online, which leads to the considerable attention of the realtime training methods of RBFNNs. Among different realtime training methods, the gradient descent (GD) based and recursive least squares (RLS) based methods have received extensive attention and been applied to various systems [7]- [10].
As one of the most classical training methods of neural networks, the GD method and its variants have achieved good performance in various offline learning scenarios, where the training sets are known a priori [11]. Different from the offline applications of neural networks, the real-time learning, as a special case of online learning, refers to the scenarios where the samples are collected one by one in chronological order. For instance, in neural network control where the RBFNN should keep sensitivity to parameter changes, the GD based update of the neural network will be carried out whenever a new sample is collected, which can also be regarded as a stochastic gradient descent (SGD) based training method with the batch size of one. Generally, the SGD based weight update law of the RBFNN in neural network control is derived from the Lyapunov stability theory based analysis and thus can guarantee the stability and closed-loop control performance of the system [12]- [14]. By contrast, the learning performance of the RBFNN in this case is very limited because the learned knowledge is closely related to the chronological distribution of the samples input to the RBFNN. Specifically, since the objective function of SGD is squared approximation error of the current sampling time, the RBFNN will mostly allocate its approximation capabilities to the latest samples while the knowledge learned from the past samples will be gradually forgotten, which is named passive knowledge forgetting phenomenon in this paper.
Due to passive knowledge forgetting phenomenon, the RBFNN needs to persistently review the past samples to realize the real-time learning (parameter convergence to the optimal value). From the perspective of adaptive control, this is the reason why the weight convergence of the RBFNN based on SGD requires the persistent excitation (PE) condition [15], [16]. When employing an RBFNN as the function approximator, what we expect is to obtain accurate approximation of the unknown function over a large-scale feature space (i.e. input space) region of the RBFNN rather than over local regions filled with the latest samples. Therefore, the training method with memory of the past samples is required to improve the real-time learning performance of the RBFNN.
As a result of its linearly parameterized form, the realtime training of the RBFNN can also adopt the RLS method. While the SGD method optimizes an objective function of the current squared approximation error, the objective function of the RLS method is the sum of squared approximation errors at each sampling time. Therefore, the RBFNN trained through RLS will allocate its approximation capabilities uniformly to each sample and keep memory to all of the past samples during the real-time learning process. Moreover, it is claimed that the RLS based training method has faster learning speed compared with the SGD based method [17]- [20]. However, as the number of samples gradually increases, the approximation capabilities of the RBFNN will be insufficient to express the excessive knowledge such that the RBFNN will be insensitive to new samples, which is also known as the data saturation phenomenon. With data saturation, the RBFNN will allocate its approximation capabilities mostly to the feature space regions where the samples are dense, resulting in limited learning performance.
To overcome the data saturation phenomenon of RLS, different forgetting mechanisms such as the forgetting factor are adopted [21]- [26]. However, while data saturation is suppressed by these forgetting mechanisms, the passive knowledge forgetting phenomenon occurs, i.e. the knowledge carried by the past samples will be gradually forgotten as the realtime learning process goes on [27]. The directional forgetting factor discussed in [19] restricts forgetting to the informationrich subspace, and suppresses the passive knowledge forgetting phenomenon to some extent. However, it is computationally complex and thus inappropriate for real-time applications. Moreover, with these forgetting factor RLS (FFRLS) methods, the learning performance of the RBFNN will still be affected by the spatial distribution of the samples. It is still possible for RBFNN to allocate its approximation capabilities too much to the regions where the samples are relatively dense, leading to poor learning of the knowledge carried by the samples within the regions where the samples are sparse.
For RLS based training methods whose objective function is the sum of the squared approximation errors of each sample, the learning performance of the RBFNN is also influenced by the distribution of samples in the feature space. To obtain the function approximation which has good generalization capability over the feature space, the samples learned by the RBFNN should be uniformly distributed over the feature space while the number of the samples should be large enough [28]. For offline learning scenarios, we can preprocess the training set such that the samples are uniformly distributed. However, for real-time applications where the samples are mostly collected in chronological order with the same sampling period, it is difficult to distribute the samples uniformly. In this case, the RBFNN cannot allocate its approximation capabilities uniformly to the feature space such that the approximation accuracy over different regions of the feature space is different.
From the aforementioned analyisis, the existing real-time training methods for RBFNNs based on SGD or RLS have the following problems to be solved: i) there is no method to suppress passive knowledge forgetting and data saturation at the same time with sufficiently little computational complexity, ii) with these methods, the RBFNN cannot allocate its approximation capabilities uniformly to the feature space, leading to the large difference in the approximation accuracy over different regions of the feature space. To address these problems, this paper proposes an RLS based real-time training method with a mechanism of memory, which is named as the selective memory recursive least squares (SMRLS) method. The design of SMRLS based training method of the RBFNN includes the following three steps: 1) The feature space of the RBFNN is normalized and then evenly discretized into a finite number of partitions, which is named as the finite-state implementation of the feature space. 2) In a real-time learning process, the memory update functions are designed according to a priori knowledge to synthesize samples within the same partition into one sample.
3) The synthesized objective function is obtained by summing up the squared approximation errors of the synthesized samples to replace the original objective function in RLS method. Similar to the ordinary RLS method, a recursive algorithm is proposed to calculate the optimal weights of the RBFNN minimizing the synthesized objective function in the real-time learning process. As an example, specific memory update functions are proposed based on two classical assumptions respectively: i) the latest samples are more reliable than the past ones [27], ii) the measurement is accurate while the function to be approxiamted is time-invariant. Based on these memory update functions, the SMRLS method is simplified such that its computational complexity is reduced to the level of the ordinary RLS method.
The discretization of the feature space and the calculation of the synthesized objective function according to memory update functions are collectively called the memory mechanism of SMRLS. Thanks to ths mechanism, SMRLS obtains the following merits over the SGD and RLS methods: 1) Through the finite-state implementation, the maximum number of the synthesized samples is set to the number of the partitions, so the data saturation phenomenon is suppressed. The knowledge of the samples within the same partition will be synthesized rather than gradually forgotten, so passive knowledge forgetting is also suppressed. This improves the learning speed and accuracy of the RBFNN. 2) Since the samples within the same partition only contribute one synthesized sample to the objective function, the RBFNN allocates its approximation capabilities uniformly to the feature space, i.e. the approximation effects over different regions of the feature space are similar. Therefore, the generalization capability of the learned knowledge is improved. The rest of this paper is organized as follows: Section II introduces preliminaries and formulates the problem. In Section III, the passive knowledge forgetting phenomenon of SGD and FFRLS methods is analyzed. Section IV provides the design process of SMRLS method in detail and analyzes its merits. Section V provides the simulation results to demonstrate the effectiveness of SMRLS. Finally, Section VI concludes the main results.

II. PROBLEM FORMULATION AND PRELIMINARIES
In Section II-A, an open-loop real-time function approximation problem is introduced. Section II-B, II-C, II-D provide a brief introduction to the RBFNN based function approximation, SGD method and RLS method, respectively.

A. Real-Time Function Approximation: An Open-Loop Case
Consider an open-loop real-time approximation task where the output of the approximator will not influence the state orbit of the system. The function to be approximated is described as follows: defined on the compact set Ω x , y ∈ R is the measurement of f (x), and v ∈ R is the measurement error including the measurement noise. The state orbit starting from the initial state x(0) is denoted by µ(x(0)) during the learning task.

Assumption 1.
To facilitate understanding of the theoretical analysis in subsequent sections, it is assumed that the measurement is accurate such that v = 0.
Assume that the output of the approximator isf (x, θ) with the parameter θ, and the following objective function is considered: Let θ * denote the optimal parameter minimizing J(θ), and we have: Let i = 1, 2, . . . , m denote the serial number of the sampling time, m denote the total number of samples, and the real-time approximation task is to constantly improve the approximation of θ * as samples (x(i), y(i)) keep being collected. The following sum of squared error (SSE) objective function is adopted to substitute J(θ) to obtain an estimation of θ * . As analyzed in [28], if the samples (x(i), y(i)) are uniformly distributed over the compact set Ω x under Assumption 1, we obatin: lim which indicates the importance of the uniform distribution of samples for the approximation. Therefore, θ * can be estimated using: if there are sufficient samples uniformly distributed over Ω x . The compact set Ω x consisting of all possible values of x is named as the feature (input) space of the approximator. In practice, the scope of Ω x is designed to be the scope over which the approximationf (x, θ) is expected to be accurate. According to the aforementioned analysis, to obtain an accurate estimation of θ * , the samples (x(i), y(i)) are required to fill the feature space Ω x , which is named as a sufficient distribution. The sufficient distribution can be satisfied using the following methods: i) for active approximation scenarios where the state orbit µ(x(0)) can be designed, µ(x(0)) is designed properly such that the feature space is filled with samples; ii) for passive approximation scenarios, although µ(x(0)) cannot be designed, we can still collect sufficient samples by long-term sampling.
Even if the sufficient distribution of samples is satisfied, for real-time approximation tasks, the approximation method (6) is still difficult to realize because: i) the direct optimization of J f (θ) is computationally complex; ii) the uniform distribution of the samples is hard to guarantee. While there is still no effective method to address ii), some real-time approximation methods have been proposed to address i) including the commonly used SGD and RLS methods [15].
Remark 1. The real-time approximation problem discussed in this section is an open-loop case where the output of the approximator will not influence the state orbit µ(x(0)). In this case, the only issue to be considered for the design of the approximation method is how to accurately approximate the function. However, when the function approximator and the original system constitute a closed-loop system, the stability of the system should also be considered. Although the closed-loop function approximation methods are more complex, they are still based on the open-loop methods. Therefore, the real-time function approximation methods discussed in this paper can also be applied to closed-loop systems with some modifications used to guarantee the stability.

B. RBFNN based Function Approximation
The RBFNN is a class of single hidden layer feedforward neural networks whose activation function is the radial basis function. The output of the RBFNN can be formulated as: where χ ∈ R l is the input vector, f N N (χ) ∈ R is the corresponding output, W = [w 1 , w 2 , . . . , w N ] T ∈ R N is the weight vector, N is the number of neurons, and Φ(χ) = [φ 1 (χ), φ 2 (χ), . . . , φ N (χ)] T ∈ R N is the regressor vector composed of radial basis functions φ i (χ), i = 1, 2, . . . , N . The radial basis function employed in this paper is the Gaussian function as follows: where c i ∈ R l is the center and σ i ∈ R is the receptive field width of the radial basis function φ i (χ).
Consider the function approximation problem in Section II-A. RBFNNs have the ability to approximate any continuous function f (x) : R n → R over the compact set Ω x ⊆ R n to arbitrary accuracy on the premise of sufficient neurons and suitably selected centers and receptive field widths [1]. Then f (x) can be approximated as: where W * ∈ R N is the optimal weight vector of the approximation, and (x) is the bounded approximation error satisfying | (x)| < * , ∀x ∈ Ω x . In this paper, the optimal value of W * ∈ R N is defined as follows: where the compact set Ω x is known as the feature space (input space) of the RBFNN. According to (4), (6), if there are sufficient samples uniformly distributed over Ω x , W * can be estimated by minimizing J f (W ) as follows: Assumption 2. It is assumed that, the RBFNN has a lattice distribution of the neuron centers, which covers the (normalized) feature space of the RBFNN.
The persistent excitation (PE) condition is an important concept in the function approximation problems based on the linearly parameterized approximators including the RBFNN and is defined as follows. [29]): The piecewise-continuous and uniformly-bounded regressor vector Φ(x) : R n → R N , is said to satisfy the PE condition, if there exist positive constants δ, α 1 , α 2 such that

Definition 1. (PE Condition
holds for ∀t 0 0, where I ∈ R N ×N is an identity matrix and α 1 , α 2 are known collectively as the PE levels.
In this paper, the amount of knowledge learned by the RBFNN is measured through the memory of samples, which is defined as follows.
where m > 0 is a small approximation error upper bound.

C. Stochastic Gradient Descent based Training for RBFNNs
To approximate W * in real-time, the following squared error objective function is considered: where k is the current sampling time. The gradient of J GD (W, k) corresponding to W is formulated as: Then W is updated in the direction of −∇J GD (W, k) as follows: where η ∈ R + is the designed learning rate.
The following theorem indicates the rationality of the weight update law (16).
Theorem 1. Assume that the sampling period is small enough. Consider the weight update law (16), if the regressor vector Φ(x) satisfies the PE condition and the learning rate η is designed properly, W will exponentially converge to a small neighborhood around W * .
Proof. Substitute (9), (15) into (16) and we obtain: whereW = W * − W is the estimation error of the optimal weight vector W * . Since the sampling period is small enough, (17) can be reformulated into the continuous form as follows: where η 0 ∈ R N ×N is the designed learning rate which is proportional to η. SinceW = W * − W , we have: Since Φ(x) is bounded and (x) is very small, system (19) can be regarded as a nominal system with a small perturbation. The nominal part of (19) is described as: According to Theorem 2.5.1 in [15], with the PE condition of Φ(x), (20) is exponentially stable. To analyze the convergence of system (19), Lemma 9.2 in [30] is considered. This nonvanishing perturbation lemma indicates that with properly designed η 0 , the small perturbation term η 0 Φ(x) (x) will not lead to large steady-state deviations ofW from zero. Therefore, under the PE condition of Φ(x), W will exponentially converge to a small neighborhood around its optimal value W * with properly designed learning rate η.

D. Recursive Least Squares based Training for RBFNNs
Different from many other neural networks, which are highly nonlinear in the parameters and require nonlinear optimization methods for learning, the RBFNN can learn with linear optimization methods such as the RLS method thanks to its linearly parameterized form.
Consider the following SSE objective function which has been mentioned in (11): Let W * LS (k) ∈ R N denote the optimal weight vector minimizing J LS (W, k), and it can be formulated as: Let R k denote the output sample vector. It is assumed that there are sufficient samples such that the matrix H T k H k ∈ R N ×N is positive definite. Therefore, the optimal weight vector W * LS (k) can be calculated as follows: According to (11), if there are sufficient samples (x(i), y(i)), i = 1, 2, . . . , m uniformly distributed over Ω x , W * can be estimated as follows: When the samples are obtained from the online sampling of the chronological data, the RLS method can be applied to estimate W * LS (k) in real-time. Let P (k) = (H T k H k ) −1 denote the covariance matrix at sampling time k, and two equivalent RLS methods are given as follows: where λ ∈ (0, 1] is the forgetting factor used to suppress the data saturation phenomenon [31]. Let P (0) = p 0 I, W (0) = W 0 denote the initial values of P (k) and W (k) respectively, where I ∈ R N ×N is an identity matrix. According to Theorem 1 in [32], the real objective function of the RLS methods (25) and (26) is formulated as: Therefore, when using (25) or (26) to update W (k), we obtain: Let λ = 1, i.e. the RLS method is applied without the forgetting factor, to make J RLS (W, k) coincide with J LS (W, k) in (21), p 0 should be chosen as large as possible.
Assumption 3. Assume that the matrix H T k H k is positive definite such that W * LS (k) exists. In this paper, the RLS method is initialized with a sufficiently large p 0 such that the weight vector W (k) based on (25) or (26) with λ = 1 satisfies: where * k > 0 is a small estimation error upper bound corresponding to the samples and RBFNN.
The following lemma shows the derivation of the ordinary RLS method which does not adopt a forgetting factor. Lemma 1. Consider the SSE objective function (21), if the following constraints: hold for ∀k ∈ N + such that H T k H k is positive definite, then W (k + 1) can be calculated using W (k) as follows: Proof. According to (23), W (k) and W (k + 1) can be formulated as: Since ∈ R k+1 , we have:

III. PASSIVE KNOWLEDGE FORGETTING AND FINITE MEMORY OF RBFNNS
The passive knowledge forgetting phenomenon of SGD and FFRLS based RBFNNs is introduced in Section III-A. By contrast, the finite memory of the ordinary RLS based RBFNN is analyzed in Section III-B.

A. Passive Knowledge Forgetting Phenomenon
The passive knowledge forgetting phenomenon can be analyzed from the perspective of objective function. For SGD training method, the objective function J GD (W, k) is shown in (14). It should be noted that the squared error term y(i) − W T Φ(x(i)) 2 only appears in J GD (W, k) at sampling time i indicating that the update of the RBFNN will not aim at reducing the error for the update at sampling time i. Therefore, even if the error y(i) − W T Φ(x(i)) 2 is reduced by the update at sampling time i, it may increase in the subsequent training process. Similarly, for FFRLS training method whose objective function J RLS (W, k) is shown in (27), the weight of the squared error term y(i) − W T Φ(x(i)) 2 , i.e. λ k−i , will gradually decrease as the real-time learning process goes on. The term "passive knowledge forgetting" is used to summarize this phenomenon in which the past samples will gradually lose the attention of the approximator as the learning process goes on. Consider the real-time approximation problem of f (x) mentioned in Section II-A. To facilitate the understanding, this passive knowledge forgetting phenomenon is presented with a special case in Fig. 1 T ∈ R 2 is the input of the RBFNN, and the state orbit of x is divided into two segments: i) an explored segment µ a and ii) a new segment µ b to be explored. Let W * a denote the optimal weight vector of the RBFNN along µ a and we have: Let Ω a denote the compact set consisting of all possible values of x along µ a . In this case, f (x) can be locally approximated as follows: It is assumed that after real-time learning along µ a using SGD or FFRLS methods, the weights of the RBFNN have converged to W * a and thus realize locally accurate approximation of f (x) along µ a . However, when the segment µ b is sampled and learned by the RBFNN, the converged weights will gradually deviate from W * a and thus lose the accurate approximation of According to Definition 2, passive knowledge forgetting will make the RBFNN lose its memory of past samples but obtain memory of latest samples. Therefore, when SGD or FFRLS methods are used to train the RBFNN, the knowledge learned by the RBFNN cannot be accumulated effectively as learning goes on, which leads to the unqualified real-time learning performance.

B. Finite Memory of the Ordinary RLS based RBFNN
For ordinary RLS method which does not adopt a forgetting factor, the squared approximation error of each sample always has the same weight in the objective function. Therefore, in real-time learning, the RBFNN will evenly allocate its approximation capabilities to each sample and tries to keep memory of as many samples as possible. The following lemma shows this memory characteristic of the ordinary RLS based RBFNN.  (25) or (26) with λ = 1 is used to train the RBFNN to approximate f (x) under Assumption 3. If the RBFNN has reasonable hyperparameter settings, then it can keep good memory of the past samples as described in Definition 2. Specifically, assume that the serial number of current sampling time is k, and then there exists a small constant * m > 0 such that the following inequation: holds for ∀i k m, where m ∈ N + is a finite integer.
Proof. Under Assumption 1, the approximation error of W * LS defined in (22) can be formulated as follows: where W * is the optimal weight vector defined in (10). Substitute (9) into (39) and we obtain: where * > 0 is the approximation error upper bound that can be arbitrarily reduced by adjusting the hypermeters of the RBFNN. Consider the difference between ordinary RLS and least squares (LS) caused by the initial value p 0 of the covariance matrix P (k) mentioned in Assumption 3, and the approximation error of the RLS based RBFNN satisfies: (41) According to (29), the following inequation is considered: where N is the number of the neurons, * k > 0 is the small estimation error upper bound caused by the initial value p 0 of P (k).
Substitute (40), (42) into (41) and we obtain: Therefore, for ∀i k m, the approximation error of the past sample (x(i), y(i)) satisfies: (44) Since * k , * can be arbitrarily reduced by designing p 0 and the hypermeters of the RBFNN while k m is a finite value, the term kN 2 * 2 k + k * 2 in (44) can be arbitrarily reduced. Let * k = kN 2 * 2 k + k * 2 and this ends the proof.
Lemma 2 indicates that the ordinary RLS based training method makes the RBFNN keep memory of every past sample and thus overcomes the passive knowledge forgetting. However, the ordinary RLS still has the following problems to be addressed, which limits the real-time learning performance of the RBFNN: 1) The approximation error upper bound in (44) is proportional to the number of samples. In practical applications, the memory of the RBFNN is finite, so when the number of the samples increases continuously, the approximation capabilities of the RBFNN will gradually become insufficient to express the excessive knowledge. Therefore, when the number of the samples is large enough, the RBFNN will almost no longer learn new samples to retain the memory of the learned samples, which is named as the data saturation phenomenon.
2) As analyzed in Section II-A, to guarantee the generalization capability of the learned knowledge, the samples included in the objective function should be uniformly distributed over the feature space, i.e. the RBFNN should allocate its finite memory uniformly to the feature space. However, with ordinary RLS, the distribution of samples is determined by the state orbit µ(x(0)) of the input vector x, resulting in the different approximation accuracy over different regions of the feature space.

IV. SELECTIVE MEMORY RECURSIVE LEAST SQUARES
To overcome the aforementioned passive knowledge forgetting and data saturation while allocating the finite memory of the RBFNN uniformly to the feature space, a real-time training method named selective memory recursive least squares (SMRLS) is proposed in this section.

A. Finite-state Implementation of the Feature Space
Although the ordinary RLS method essentially minimizes the objective function (27) with λ = 1 using a recursive method, Lemma 1 indicates that the derivation of RLS actually considers optimization of the SSE objective function (21). Therefore, to obtain a training method uniformly allocating the finite memory of the RBFNN to its feature space, the following SSE objective function is considered: where k is the serial number of the current sampling time, M (k) is the number of samples at sampling time k, and the sampled input x(j), j = 1, 2, . . . , M (k) is uniformly distributed over the feature space Ω x for ∀k ∈ N + . The most prominent characteristic of J LS (W, M (k)) is that its samples are always uniformly distributed over the feature space during the real-time learning process.
To obtain an objective function in the form of (45), the following steps, which are collectively named as the finitestate implementation of the feature space, are adopted: 1) The feature space Ω x of the RBFNN is normalized into a certain scope Ωx. For example, Fig. 2 shows a twodimensional case of the normalization into the scope of 1]. In this case, the original input of the ∈ Ω x is normalized as follows: T ∈ Ωx is the normalized input. The reason for the normalization is to simplify the neuron center distribution of the RBFNN such that f (x) along different state orbits can be approximated using the same neuron distribution satisfying Assumption 2 as follows: whereW * is the optimal weight vector defined as: 2) The normalized feature space Ωx is evenly divided into N P partitions which are mutually disjoint as shown in Fig. 2. Let Ω j x , j = 1, 2, . . . N P denote the N P partitions, and the samples within each partition Ω j x are synthesized into one sample (γ j (k), ϕ j (k)) in real-time. Assume that there are M P (k) partitions which have been sampled before sampling time k and the number of partitions N P is large enough, if the synthesized samples satisfy γ j (k) ∈ Ω j x , j = 1, 2, . . . , M P (k), the uniform distribution of the synthesized samples is guaranteed. To facilitate the analysis, in a certain real-time learning process, the serial number j of a partition is assumed to be the order in which this partition is first sampled. For example, the partition in which the first sample (x(1), y(1)) is collected is denoted by Ω 1 x . With the finite-state implementation, the number of the synthesized samples at sampling time k is equal to the number of the sampled partitions M P (k). Therefore, at sampling time k + 1, the number of the synthesized samples M P (k + 1) satisfies: where Ω j x ⊂ Ωx is the union set of the partitions which have been sampled at least once before sampling time k. Therefore, the SSE objective function of the synthesized samples can be formulated as follows: which is named as the synthesized objective function and coincides with the form of (45).
Since (γ j (k), ϕ j (k)) represents the synthesis of the samples within the partition Ω j x at sampling time k, its value will be updated only when the samples within Ω j x increase. Therefore, a partition judgment function is defined as follows: The partition judgment function G(j,x(k)) is used to determine in which partition the current sample (x(k), y(k)) locates and thus to determine which synthesized sample should be updated. Based on the partition judgment function, the update laws of γ j (k) and ϕ j (k) are proposed as follows: γ j (k + 1) = γ j (k) + G(j,x(k + 1))S γ (γ j (k),x(k + 1)) ϕ j (k + 1) = ϕ j (k) + G(j,x(k + 1))S ϕ (ϕ j (k), y(k + 1)) (52) where S γ (·) and S ϕ (·) are memory update functions representing the update mechanism of the synthesized samples. Specially, if the new knowledge (x(k+1), y(k+1)) is believed to be more reliable than the old knowledge (γ j (k), ϕ j (k)), the memory update functions S γ (·), S ϕ (·) can be designed as: Remark 2. It should be noted that the implementation of the memory update law (52) requires the initialization of the synthesized sample (γ j (k), ϕ j (k)). In this paper, the initial value of (γ j (k), ϕ j (k)) is set to the first sample within the partition Ω j x during the online learning process. Let k j k denote the sampling time when partition Ω j x is first sampled, and the initialization of (γ j (k), ϕ j (k)) is described as follows: where j = 1, 2, . . . , M P (k) represent the serial numbers of the partitions which have been sampled at least once before sampling time k.
Through the normalization and partitioning operations, the finite memory of the RBFNN is uniformly allocated to the normalized input space Ωx and the maximum quantity of the synthesized samples is limited to the total number of the partitions, i.e. N P . With the finite-state implementation, the aforementioned uniform sample distribution and data saturation problems of the ordinary RLS method are addressed.

B. Recursive Optimization for Synthesized Objective Function
To calculate the weight vector W (k) minimizing J S (W, M P (k)) in real-time with little computational complexity, a method named selective memory recursive least squares (SMRLS) is proposed. The derivation of SMRLS is shown in detail as follows.
Assumption 5. It is assumed that the sample (x(k + 1), y(k + 1)) is collected from the partition Ω ā x at sampling time k + 1. Consider the objective function (50) and the memory update functions (52). Different from the SSE objective function of the ordinary LS method in (21) whose samples continuously increase at each sampling time, the objective function (50) of SMRLS is about the synthesized samples (γ j (k), ϕ j (k)), j = 1, 2, . . . , M P (k) whose number does not necessarily increase as the learning process goes on. At the sampling time k + 1, the behaviors of J S (W, M P (k)) can be classified into two categories: addition or replacement. Under Assumption 5, these two behaviors can be explained as follows: 1) The addition behavior means that if the sample (x(k + 1), y(k + 1)) is the first sample within the partition Ω ā x , a new synthesized sample (γ a (k + 1), ϕ a (k + 1)) is generated and a new term ϕ a (k + 1) − W T Φ(γ a (k + 1)) 2 will be added to the synthesized objective function (50). In this case, a = M P (k + 1) = M P (k) + 1.
2) The replacement behavior means that if the sample (x(k + 1), y(k + 1)) is not the first sample within the partition Ω ā x , the synthesized sample (γ a (k), ϕ a (k)) will be updated to (γ a (k + 1), ϕ a (k + 1)) according to (52) and the squared error term ϕ a (k) − W T Φ(γ a (k)) 2 in (50) will also be replaced by a new term ϕ a (k + 1) − W T Φ(γ a (k + 1)) 2 . In this case, a M P (k + 1) = M P (k). To distinguish these two behaviors, a function of the latest sampled input vectorx(k + 1) is defined as follows: where Ω M P (k) is the union set of the partitions which have been sampled at least once before sampling time k.
To obtain the update law of SMRLS, the update from sampling time k to k + 1 is analyzed. According to (50), the SMRLS based weight vector W (k) and W (k+1) are assumed to satisfy the following constraints: Based on the aforementioned addition and replacement behaviors of J S (W, M P (k)), the update from W (k) to W (k + 1) can be classified into two cases as follows: 1) Ifx(k + 1) / ∈ Ω M P (k) , which means that the sample (x(k + 1), y(k + 1)) is the first sample within the partition Ω ā x , the aforementioned addition behavior occurs. Therefore, we have a = M P (k + 1) = M P (k) + 1, and the synthesized objective function J S (W, M P (k + 1)) can be described as: Consider (56), (57), and the constraints of the update can be summarized as: which coincides with the constraints of the ordinary RLS as analyzed in Lemma 1. Therefore, ifx(k + 1) / ∈ Ω M P (k) , the update from W (k) to W (k + 1) can adopt the ordinary RLS based update law (31), (32) and be formulated as follows: S (k + 1) = P −1 S (k) + Φ(γa(k + 1))Φ T (γa(k + 1)) Φ(γa(k + 1)) = [φ1(γa(k + 1)), . . . , φN (γa(k + 1))] T (59) where P S (k) is the corresponding covariance matrix. According to the analysis in Section II-D, the initial value of P S (k) is set to p 0 I, where I ∈ R N ×N is an identity matrix and p 0 > 0 is as large as possible.
According to Remark 2, the value of (γ a (k + 1), ϕ a (k + 1)) is set to (x(k + 1), y(k + 1)), and thus (59) becomes: (60) 2) Ifx(k + 1) ∈ Ω M P (k) , which means that (x(k + 1), y(k + 1)) is not the first sample within the partition Ω ā x , the aforementioned replacement behavior occurs. Therefore, we have a M P (k + 1) = M P (k), and the synthesized objective function J S (W, M P (k + 1)) can be described as: To obtain the recursive update law of W (k), a virtual synthesized objective function, which involves the synthesized samples of the partitions apart from Ω ā x , is proposed as follows: Since J V S (k + 1) = J V S (k), substitute (62) into (61), and we have: Similar to J S (W, M P (k + 1)), the synthesized objective function J S (W, M P (k)) can be described as follows: denote the virtual weight vector minimizing the objective function J V S (k) as follows: Consider (56), (63), (64), (65), and we have the following update constraints: J S (W, M P (k + 1)) = J V S (k + 1) + ∆J S (k + 1).

Remark 3.
It should be noted that the implementation of (70) requires storage of the past synthesized samples of each partition. At sampling time k + 1, each partition Ω j x , j = 1, 2, . . . , M P (k + 1) should hold the memory of (γ j (k), ϕ j (k)) to realize the recursive calculation. Assume that the number of partitions N P is large enough such that Φ(γ a (k + 1)) − Φ(γ a (k)) is very small and Φ(γ a (k)) can be approximated by Φ(γ a (k + 1)). Substitute Φ(γ a (k)) = Φ(γ a (k + 1)) into (70), and we obtain: (72) which only requires the storage of ϕ j (k). To reduce the computational complexity, (72) can be used to replace (70) in the recursive optimization process. The rationality of this simplification is based on the premise that the number of partitions N P is large enough, so when N P is small, this simplification will reduce the learning accuracy of SMRLS.

Remark 4.
Since the maximum number of the samples in SMRLS is limited to N P , i.e. the number of partitions, the accuracy of SMRLS is related to the partitioning operation. When N P is too small, the accuracy of SMRLS will be affected while too many partitions will increase the computational complexity. Therefore, N P should be designed properly according to specific learning tasks and partitioning methods other than evenly partitioning deserve further exploration.

C. SMRLS Based on Common Assumptions
The implementation of SMRLS requires the specific design of the memory update functions in (52) based on a priori knowledge of the system. Different memory update functions will lead to very different characteristics of the SMRLS based methods. For real-time function approximation scenarios, there is a classical assumption that the latest samples are more reliable than the past ones, which is also the basis of the SGD and FFRLS methods [27]. Interestingly, SMRLS can also be simplified with this assumption and the derivation of the simplified method is shown in detail as follows.
Assumption 6. In the SMRLS based real-time function approximation process, the latest sample within the partition Ω j x , j = 1, 2, . . . , M P (k) can represent the unknown function more accurately compared with the past samples within the same partition.
(76) Similar to the analysis of Remark 3, when N P is large enough, SMRLS-R can also be simplified as follows: (77) which reduces the computational complexity of SMRLS-R.
In many practical systems, the unknown functions to be approximated are time-invariant, and thus the following assumption is also considered.
Assumption 7. The function to be approximated by the RBFNN is time-invariant and the measurement is accurate.
Under Assumption 7, each sample within the same partition is regarded to have the same effect to characterize the unknown function. Therefore, the first sample within the partition Ω j x , j = 1, 2, . . . , N P is taken as the memory of Ω j x and not updated in the subsequent learning process. In this case, the memory of each partition is fixed after the learning of the first sample within the partition, and thus the method is named as fixed SMRLS (SMRSL-F). Consider (71), and we only need to stop the update when ρ(x(k+1)) = 1 to obtain the description of SMRLS-F: ∈ Ω M P (k) : 2) if ρ(x(k + 1)) = 1, i.e.x(k + 1) ∈ Ω M P (k) : (78) In this case, the memory update functions in (52) are regarded as S γ (γ j (k),x(k+1)) = S ϕ (ϕ j (k), y(k+1)) = 0 after the initializatin of (γ j (k), ϕ j (k)). Compared with the ordinary RLS method, SMRLS-F obtains lower computational complexity while its accuracy is influenced by the number of partitions according to Remark 4.

D. Merits of Selective Memory Recursive Least Squares
Compared with the classical RLS or SGD based training methods, SMRLS has the following merits: 1) Selective Memory of the RBFNN: Since the derivation of SMRLS is based on the ordinary RLS method, the SMRLS based RBFNN also has the finite memory mentioned in Lemma 2 such that all of the knowledge (samples) will be synthesized and learned by the RBFNN during the real-time learning process. In this way, with properly designed memory update functions, the passive knowledge forgetting phenomenon is effectively suppressed.
In addition, the finite-state implementation of the feature space also suppresses the data saturation problem of the ordinary RLS method. According to (71), if ρ(x(k + 1)) = 0, the eigenvalues of P −1 S (k + 1) satisfies: where λ i (P −1 S (k)), i = 1, 2, . . . , N represent the eigenvalues of P −1 S (k) and tr(·) is the trace operator. Since P −1 S (k) is positive definite, i.e. λ i (P −1 F (k)) > 0, the eigenvalues of P −1 S (k) will gradually increase as more partitions are explored and thus the eigenvalues of P S (k) will gradually decrease. The term P S (k+1)Φ(x(k+1)) in case 1) of (71), which represents the learning rate, satisfies the following constraint: PS(k + 1)Φ(x(k + 1)) 2 λ 2 max (PS(k + 1)) Φ(x(k + 1)) 2 N λ 2 max (PS(k + 1)) = N/λ 2 min (P −1 S (k + 1)) (80) where N is the number of the neurons, λ max (·) and λ min (·) represent the maximum and minimum of the eigenvalues, respectively. According to (71), the learning rate P S (k + 1)Φ(x(k + 1)) will gradually decrease as more partitions are explored, which means that the finite memory of the RBFNN is being consumed. With the ordinary RLS method, this memory consumption happens whenever a new sample is learned by the RBFNN and thus the weights W (k) will nearly stop updating when the finite memory is depleted, i.e. the data saturation occurs. By contrast, with the SMRLS method, the maximum number of the synthesized samples is set to N P , so the finite memory of the RBFNN will not be depleted during the real-time learning process. Therefore, SMRLS also effectively suppresses the data saturation phenomenon.
Since both the passive knowledge forgetting and the data saturation are suppressed through the selective memory, the real-time learning speed and accuracy of the RBFNN are improved.

2) Uniformly Allocated Approximation Capabilities:
As analyzed in Section II-A, the uniform distribution of the samples over the feature space is essential for the generalization capability of the learned knowledge. With SMRLS, if the number of partitions N P is large enough, the synthesized samples (γ j (k), ϕ j (k)), j = 1, 2, j = 1, 2, . . . , M P (k) can be regarded to be uniformly distributed over the feature space. Therefore, the SMRLS based RBFNN allocates its finite approximation capabilities uniformly to the feature space and thus obtains good generalization capability.

V. SIMULATION STUDIES
In this section, the RBFNN is employed to approximate an unknown function of a single inverted pendulum-cart system, whose state space description is formulated as follows: .
The state variables x 1 and x 2 represent the angular position and velocity of the pendulum, respectively, g = 9.8m/s 2 is the gravitational acceleration, m c = 0.1kg is the weight of the cart, m = 0.02kg is the weight of the pendulum, and l = 0.2m is half the length of the pendulum. The RBFNN is used to approximate f I (x) in real-time.
In this section, the following real-time learning methods are compared: SMRLS-R (76), ordinary RLS (25)  It should be noted that this paper focuses on the realtime learning problem of the RBFNN rather than the control. Therefore, the RBFNN is not involved in the closed-loop system and it is assumed that the output of the unknown function f I (x) to be approximated can be measured directly as follows: The objective of the real-time learning task is to obtain an accurate approximation of f I (x) over the normalized feature space Ωx using the increasing samples (x(k), y(k)).

A. Simulation Studies for Data Saturation
To test the suppression effect of different learning methods to data saturation, the following two real-time learning scenarios are considered.

(A) Abrupt parameter perturbation
Consider the following state orbit of x 1 (t) and corresponding measurement: where the value of l in (82) is changed abruptly as below: The training of RBFNNs is carried out over [0, 100s]. Fig. 3 and Fig. 4 demonstrate the convergence process of the approximation errors and weights, respectively. In the first 50s of the training, the RBFNNs have approximated the unknown function f I (x) with the parameter l = 0.2, and thus when the parameter abruptly changes, the data saturation phenomenon greatly reduces the convergence speed of the ordinary RLS method. Compared with SGD and FFRLS methods which can suppress data saturation, SMRLS-R achieves faster learning speed to the changed function. According to the derivation of SMRLS in Section IV, the SMRLS-R based weight convergence can be completed within only one period of the state orbit (2πs for case (A)) if the sampling period is small enough such that all of the partitions along the state orbit can be sampled at least once in a single period.
To test the accuracy of the approximation, the weights at 100s are recorded as the learned knowledge to approximate f I (x) with the latest parameter l = 0. 3. Let e f = f I (x) − W T Φ(x) denote the approximation error. As shown in Fig. 5, while all the four methods achieve the basic approximation of f I (x), the accuracy of FFRLS and SMRLS-F methods is obviously better than the other two methods. The reason why SMRLS-R is less accurate than FFRLS lies in the limitation of the maximum number of samples, i.e. N P , mentioned in Remark 4. The increase of N P will improve the approximation accuracy but also increases the computational complexity and the requirement for sampling period.

(B) Decreasing measurement error
Consider the following state orbit of x 1 (t) over [0, 100s] and corresponding measurement: where v(k) = e −0.2k sin k is the decreasing measurement error. Case (B) represents the scenarios where the output of the unknown function cannot be measured directly and requires extra estimation with the converging measurement error. For example, in a class of control systems, the output of the unknown function to be approximated can only be obtained using an observer and the objective of the RBFNN is to approximate the function with the measurements. Fig. 6 shows the convergence of the approximation error in the real-time training process, indicating that SMRLS-R still has the fastest convergence speed. The approximation effect using the learned knowledge is shown in Fig. 7, where FFRLS has the highest learning accuracy and the accuracy of SMRLS-R is higher than the ordinary RLS method.
According to the results in case (A) and (B), SMRLS-R suppresses the data saturation phenomenon effectively, and thus has faster convergence speed for time-varying learning scenarios compared with the other classical real-time training methods. However, the accuracy of its approximation is determined by the partitioning operation.

(C) Non-repetitive state orbit
To test the suppression effect of SMRLS to passive knowledge forgetting, the following non-repetitive state orbit is considered: where the duration of the learning task is set to 100s.
Since case (C) coincides with Assumption 7, SMRLS-F is also compared with the other methods in this case. Fig. 8 shows that FFRLS has the highest real-time tracking accuracy of f I (x) along the non-repetitive state orbit while SMRLS-R, ordinary RLS and SMRLS-F have similar tracking accuracy. However, high tracking accuracy does not mean high accuracy of the learned knowledge. As shown in Fig. 9, the weights of FFRLS maintain a relatively large updating speed over the training process indicating that its high-accuracy tracking of f I (x) may be at the cost of forgetting the past samples. The weight vectors at 100s are recorded as the learned knowledge and Fig. 10 demonstrates the learning accuracy of the methods. Consistent with the analysis in Section III, the passive knowledge forgetting of FFRLS and SGD improves their real-time tracking accuracy of f I (x) while reducing the accuracy of the learned knowledge because the past samples are gradually forgotten. By contrast, the SMRLS and ordinary RLS methods sufficiently utilize both the past and the latest samples such that the accuracy of the learned knowledge is guaranteed. Moreover, the results also indicate that SMRLS-F has very similar learning performance to SMRLS-R in case (C) and thus becomes a computationally efficient alternative of SMRLS-R under Assumption 7.
C. Simulation Studies for Approximation Capability Allocation (D) Non-uniformly distributed state orbit One of the most important merits of SMRLS is its uniform allocation of the approximation capabilities such that the RBFNN will pay average attention to its approximation effect in different regions of the feature space and thus obtains good generalization capability. To test this merit, the random Non-Uniform Rational B-Splines (NURBS) state orbit µ(x(0)) in Fig. 11(a) is adopted and the value of f I (x) along this orbit is shown in Fig. 12.
Since the learning of SGD is less effective than the RLS based methods and SMRLS-F has very similar learning effect to SMRLS-R, only SMRLS-R, ordinary RLS and FFRLS methods are compared in case (D). Fig. 11(b)-(d) and Fig. 13 demonstrate the update of the weight vectors and approximation errors, respectively. In the learning process, FFRLS has higher real-time tracking accuracy to f I (x) than SMRLS-R and ordinary RLS methods. In real-time approximation tasks including case (D), the stability of the weights reflects the generalization capability of the learned knowledge. Specifically, if the weights hardly update over a time interval, it means that this set of weights can approximate f I (x) accurately during this interval. Therefore, the dramatical update of the weights shown in Fig. 11(d) indicates the weak generalization capability of FFRLS. Fig. 13 demonstrates the accuracy of the learned knowledge along the same random NURBS state orbit. The four subplots (i) − (iv) of Fig. 13(b) are corresponding to the four segments in red of Fig. 11(a). While ordinary RLS mostly allocates the approximation capabilities of RBFNN to the regions where the samples are dense (i.e. the segments in blue of Fig. 11(a)) and FFRLS mostly allocates the approximation capabilities to the latest samples, SMRLS-R allocates the approximation capabilities uniformly to the feature space. Therefore, the FFRLS based approximation is more accurate over the regions filled with the latest samples than the other regions, the ordinary RLS approximates f I (x) more accurately over the regions where the samples are dense, and SMRLS-R has relatively uniform approximation accuracy along the state orbit. Particularly, SMRLS-R effectively improves the approximation accuracy over the regions (i) − (iv) in Fig.  11(a) where the samples are sparse compared with FFRLS and ordinary RLS methods.
The aforementioned simulation results indicate that the SMRLS method can suppress data saturation and passive knowledge forgetting at the same time and thus improves the real-time learning performance of the RBFNN. Moreover, through SMRLS, the finite approximation capabilities of the RBFNN are uniformly allocated to the feature sapce such that the learned knowledge has good generalization capability.

VI. CONCLUSION
In this paper, the real-time learning problem of RBFNNs is studied and a phenomenon named passive knowledge forgetting is analyzed to explain the limited learning performance of RBFNNs caused by the forgetting mechanism in SGD and RLS based training methods. To suppress both passive knowledge forgetting and data saturation, the SMRLS method is proposed by discretizing the feature space and optimizing a synthesized objective function with the recursive method. Through SMRLS, the RBFNN can allocate its finite approximation capabilities uniformly to its feature space and thus obtains improved real-time learning performance.
In addition to the RBFNN, SMRLS is also applicable for real-time training of other single hidden layer feedforward neural networks. Moreover, as a result of its good learning performance and acceptable computational complexity, SMRLS has considerable potentials in the field of adaptive control, system identification and machine learning.


extfieldsofstudy:
Computer Science
Physics

text:
LightAMR format standard and lossless compression algorithms for adaptive mesh refinement grids: RAMSES use case

The evolution of parallel I/O library as well as new concepts such as 'in transit' and 'in situ' visualization and analysis have been identified as key technologies to circumvent I/O bottleneck in pre-exascale applications. Nevertheless, data structure and data format can also be improved for both reducing I/O volume and improving data interoperability between data producer and data consumer. In this paper, we propose a very lightweight and purpose-specific post-processing data model for AMR meshes, called lightAMR. Based on this data model, we introduce a tree pruning algorithm that removes data redundancy from a fully threaded AMR octree. In addition, we present two lossless compression algorithms, one for the AMR grid structure description and one for AMR double/single precision physical quantity scalar fields. Then we present performance benchmarks on RAMSES simulation datasets of this new lightAMR data model and the pruning and compression algorithms. We show that our pruning algorithm can reduce the total number of cells from RAMSES AMR datasets by 10-40% without loss of information. Finally, we show that the RAMSES AMR grid structure can be compacted by ~ 3 orders of magnitude and the float scalar fields can be compressed by a factor ~ 1.2 for double precision and ~ 1.3 - 1.5 in single precision with a compression speed of ~ 1 GB/s.

Introduction
RAMSES is a massively parallel hydrodynamical code for self-gravitating magnetized flows widely used in the computational astrophysics community [1], to solve the evolution of dark matter, stellar populations, and gas via gravity, hydrodynamics, radiative transfer, and non-equilibrium radiative cooling/heating. For hydrodynamics, it uses the HLLC Riemann solver [2] and the MinMod slope limiter to construct gas variables at cell interfaces from their cell-centered values. The dynamics of collision-less dark matter and star particles are evolved with a multi-grid particle-mesh solver and cloud-in-cell interpolation [3].
It implements an adaptative mesh refinement (AMR) technique to optimize the memory/precision ratio in numerical simulations, [4], [5]. Its data structure is a cell-based fully threaded octree [6], which means that each subdomain describes a full octree of the entire simulation domain up to its topmost root grid which has the size of the whole simulation domain. In this distributed AMR tree approach each subdomain defines its own local grid refinement. In most use cases, the domain decomposition is based on a 1D Hilbert curve to split the calculation between MPI processes. It was first designed for cosmological and galaxy numerical simulation, and was used for large simulation [7], [8], [9]. However, RAMSES faces I/O scalability bottlenecks when used over a few thousand of MPI processes, Email addresses: loic.strafella@cea.fr (L.Strafella), damien.chapon@cea.fr (D.Chapon) partly due to its "multiple file per process" I/O strategy. A typical cosmological simulation using 10,000 MPI processes will produce at least 40,000 files (and filesystem inodes) per output with small file size from dozens of megabytes to a few hundred megabytes, which is far from optimal on modern Lustre filesystems and generally rapidly reach the user's limits on supercomputer. What is more, since RAMSES only implements one output data format, for both checkpointing and post-processing purposes, it is not optimized for both cases. Therefore, a simulation like Extreme-Horizon [9] produced 3.6 terabytes of data with more than 50,000 files per snapshot (AMR files, hydrodynamical quantity files, particles files, etc). Considering that several snapshots are required for scientific analysis that kind of simulation becomes rapidly challenging in terms of I/O, data management and post-processing.
In a previous work, the integration of the the Hercule library [10] for parallel I/O and data management in RAMSES was the first step to improve the I/O performance, scalability and data management [11]. In terms of data management, the integration of a parallel I/O library was also the occasion to add a post-processing specific data-flow to RAMSES. This was an important step which permitted to optimize significantly the checkpoint dump since the checkpoint's raw data will no longer be used by post-processing tools for simulation data analysis. With the original version of RAMSES, sometimes the size of the simulation is calculated not based on CPU performance, memory consumption, or pure scalability of the code but simply based on the user's disk storage available capacity. We expect this problem to occur more and more often in the coming years. The usually adopted solution is either to reduce the size of the simulation, the precision of the simulation by tuning AMR parameters or reducing the number of produced snapshots. The question was then is there a standard and compact data model that can be used to describe the AMR data from RAMSES new post-processing data-flow ?
Hercule library's post-processing database strategy is to propose an XML dictionary containing the description of standardized data model for different kind of meshing strategies and objects. Therefore, once a standardized data model is described in the dictionary, every tool using Hercule will be able to "understand" Hercule post-processing data. As a reminder, we show in  Until a few years ago, AMR data did lack of a standardized format for post-processing and visualization purposes partly due to the different flavor of AMR strategies: block-based / patchedbased [12], [13], [14], cell-based (RAMSES), etc. Thus AMR data still required some data transformation step (e. g. downcasting to unstructured grids, resampling to fixed resolution cartesian grids) to be manipulated by general purpose post-processing or visualization tools (e. g. Paraview 1 , VisIt 2 ). An AMR mesh described using unstructured grids can easily lead to memory issues and produce high volume data files, thus each AMR code uses a custom data model to output the data. In the case of RAMSES, the AMR internal 32 bits integer linked-lists are dumped level by level. Code-specific data-processing libraries were then implemented as solutions to read the RAMSES datafile format and handle its AMR grid structures (e. g. PyMSES 3 , Osiris 4 ) but these libraries needed to be updated to follow the RAMSES data format evolutions [15]. A standardized and self-consistent data model aims at reducing the cost of postprocessing tools maintenance and encourage sharing development efforts.
In this work, we introduce a new lightweight AMR data model that we applied on RAMSES AMR data. By using this data model, we show that data volume of RAMSES outputs can be significantly reduced, without loss of information . In addition, we present two lossless data compression algorithms that can be built upon this new data model standard to further reduce 1 Paraview.org 2 VisIt home page 3 http://irfu.cea.fr/Projets/PYMSES 4 https://www.nbi.dk/ nvaytet/osiris/osiris.html both the size of the AMR grid refinement description and the size of associated physical quantity scalar field data (float 32/64 bits). In the results section, we present lightAMR data volume reduction as well as tree pruning and compression benchmarks on published RAMSES simulation datasets.
2. Scientific analysis-specific AMR light tree data model

Motivation
The RAMSES checkpoint/restart data format was designed, as in many other simulation codes used in astrophysics, to hold all the information necessary for the code to restart a simulation run and resume the dynamical evolution of the numerical model. Since our goal is to design a new output data format for the specific purpose of scientific analysis and post-processing, we are free to select from all the information available in memory at runtime what is strictly necessary to answer one's postprocessing need.
To fit one's post-processing needs, a user can choose to export a subset of physical quantities (e. g. store only the density and thermal pressure scalar fields defined on the AMR grid, and discard the gravitational potential/acceleration AMR fields, the velocity AMR vector field and the Lagrangian particle information), hence reducing the overall dataset volume on disk. In RAMSES, this user-defined data field subset selection can be configured at runtime using the input parameter namelist, without the need to recompile the code. What is more, for some data post-processing use cases, down-casting these quantities from double to single precision could even prove sufficient, reducing further the data volume by a factor of 2. Finally, since this new data format is not meant to be read by RAMSES itself upon restart but only by data post-processing tools, we can decide to adopt a new standardized AMR grid data structure and a custom encoding to be able to answer any type of post-processing software requirement in an optimized and standardized way.

Implementation
The Hercule I/O library HDep standardized data model for AMR trees aims to offer a light and self-consistent way to describe an AMR grid. In this work, we adopt a similar approach as the one proposed in the HDep format, the details of our approach is described here. We further reference this standardized data model for AMR grids as lightAMR. Since that kind of mesh can be expressed in the form of a tree, it uses a boolean description of the AMR tree from top to bottom, level by level, as shown on figure 1. With this breadth first traversal description of an AMR tree, two main arrays are used: one for the grid refinement description and one for the cell ownership mask. The first array describes the refinement state of a cell: 1 for a refined cell and 0 for a leaf (unrefined) cell. All child cells from a refined one must be described, as shown in 1. This implies that the number of cells at a level l must be equal to the number of refined cells at level l − 1 times the refinement factor to the power of the dimension. In our 2D example, the refinement factor being 2 by dimension, the level 2 must contain 3 × 2 2 boolean values since there are 3 refined cells at level 1. The value in the ownership mask array describes for each cell if it does not belong to the current domain (according to the domain decomposition policy): 0 if it belongs to the current domain, otherwise 1. Figure 1: Description of a 2D AMR grid (with a refinement factor of 2) dispatched on two MPI processes according to the lightAMR standardized data model (above). Grid refinement, ownership mask and scalar data field arrays for the domain #1 (below).
In order to be self-consistent, a few other (lightweight) metainformation such as the size of the simulation box, the physical quantity units, physical model details or used numerical schemes need to be stored, using an explicit semantic (string key-value pairs). In addition, some useful meta-information like the number of cells per level, the total number of grids, the number of levels, if somewhat redundant with the AMR grid description, can be conveniently stored to enhance the postprocessing tools performance. AMR-bound physical quantities are described in scalar field data arrays that follow the same ordering as the grid refinement array, hence providing a straightforward mapping (which facilitates the LOD approach). Since RAMSES describes the AMR tree from the root node which is the simulation box itself on each domain and stores scalar field values even for refined (coarse) cells, the construction of these arrays is straightforward when converting in the lightAMR format. Of course, this lightAMR tree data format can be easily extended to be compatible with other numerical codes that do not share RAMSES fully-threaded tree data structure, refined cell (coarse) scalar field values only requiring to be computed on-the-fly upon data export.
Collision-less particles are treated in a Lagrangian way within RAMSES and are linked to the AMR grid. The particles and their associated properties are simply added to our new data output format in the form of one dimensional arrays. Particle positions for example are in the form x 1 , y 1 , z 1 , ..., x n , y n , z n . We chose not to associate particles to their linked AMR cells and consequently not overload the data format. The particle-cell binding is left to the post-processing tools to be performed onthe-fly if needed.
For simulation codes that do not use the fully-threaded approach and naturally describe collections of trees at a defined minimum level l 0 on simulation domains, the LightAMR data model can be trivially extended. One just needs to provide an additional AMR root description array to detail the I-J-K logical position indices of the root coarse cells of those trees in a structured grid at level l 0 . While providing an array containing the logical position in (i 1 , j 1 , k 1 , i 2 , j 2 , k 2 , ...) format would be the natural approach, one may choose to provide an array of Morton indices of logical position at the level l 0 . Indeed, using the Morton curve, 3D/2D logical coordinates can be mapped on 1D Morton curve by bits interleaving, hence reducing the size of the array by a factor 3 (resp. 2) in 3D (resp. 2D) but with a limitation on the maximum value of l 0 32 bits Morton indices: l 0,max = 10 in 3D (l 0,max = 16 in 2D). Within this extension, compared to their fully-threaded tree version, the grid refinement, ownership mask and scalar field data arrays should be trimmed up to level l 0 as a consequence.

Tree pruning to remove redundancy
Due to the fully threaded octree approach in RAMSES, some AMR cells may be described several times in the trees of different domains, leading to data redundancy. Even if those multiple AMR cell descriptions are required by a numerical solver like the Poisson equation multi-grid solver [3], describing those cells only once is sufficient for post-processing and visualization purposes. Therefore, a tree pruning algorithm can be implemented to reduce the memory footprint of our file format.
Once the tree is built according to lightAMR format described above, most of the data redundancy in the AMR dataset can be removed. We developed a tree pruning algorithm, propagating from the bottom of the tree to its top, in order to remove the redundant cells. The unnecessary cells are removed by changing their refinement value (C.R.V on figure 2) and removing their described children cells on the next level. The C.R.V. is done only for a refined cell if all its children belong to another domain (mask = 1), as shown on figure 2. Then those changes are propagated to the associated arrays: the ownership mask and scalar field data arrays. Therefore, removing those refined cells and its children is the first step to reduce the total volume of data before applying any compression algorithm. Depending on the simulation, this step can prune a significant part of any RAMSES simulation dataset, see results on astrophysical simulation datasets in section 5.1. Prior to the tree pruning step, a temporary mask array needs to be adapted from the original ownership mask array to properly compute the C.R.V. This slight adaptation insures that any refined cell is not flagged as masked if at least one of its children is not masked. This temporary mask is only used during the tree pruning step .
At the end of this step, the AMR grid is properly described with one or several trees and most of the data redundancy was removed. However, for a fully-threaded octree the redundancy cannot be entirely removed unless the tree is split at the coarse level of refinement (a.k.a. levelmin in RAMSES) and domain boundaries set between sibling cells at that coarse level (and not deeper in the tree). The next stage in the lightAMR formatting is to apply data compression on both the grid refinement/ownership mask arrays and the simple/double precision associated scalar field data arrays.

AMR data compression
To overcome the I/O bottleneck, in addition to the integration of parallel I/O [11] and a data reduction strategy, data compression can reduce the volume of data to transfer or store on disk/object store and further improve the overall I/O performance. Combining those different strategies can have a significant impact on the I/O times, data volume storage and data transfer rates, whether for in-transit or in-situ data post-processing. In this section, we present two lossless compression schemes specific to octree AMR meshes that we developed for the ligh-tAMR data model; one for the grid description boolean valued arrays and one for the scalar field float (32/64 bits) data arrays.

LightAMR grid boolean arrays lossless compression
Once the lightAMR structure has been described with boolean valued arrays, we developed a deterministic, lossless and performant algorithm to compress the information contained in the two main arrays (grid refinement and ownership mask arrays) based on a custom run-length encoding which is well known compression technique [16] . The first step of the algorithm is to consider the boolean array as an ordered list of packs of identical bit values, either successive zeroes or successive ones, and count the lengths of these continuous packs. The second step is to digitize these ordered pack sizes with a special encoding and store the result in a single byte integer array. Out of the 256 integer values available, our encoding only need values between 1 and 62. To encode pack sizes, we use a base-52 decomposition and we also record the number of digits (only if greater than 1) in that decomposition. The uint8 values ranging from 2 to 7 are used to store the number of bytes used for the encoding of a each pack size, and the 52 uint8 values between 11 and 62 (values 9 and 10 are used for specific purposes, see next paragraph) to store the digits of the base-52 decomposition (with an additional shift of 11). For example, to encode a pack size of 2904 which decompose in 1 · 52 2 + 3 · 52 1 + 44 · 52 0 , the number of digit (3) is stored, followed by the digits [1,3,44], encoded as [12,14,55] ([1 + 11, 3 + 11, 44 + 11]), resulting in the following 4 bytes [3,12,14,55] instead of the 262 bytes required to store 2904 zeroes or ones. With that encoding, the greatest pack size that can be encoded is 52 7 − 1 ∼ 10 12 continuous 0 or 1, which is more than enough, the maximum number of cells per domain in a typical RAMSES simulation run never exceeding a 10 7 due to memory limitations. Since the encoded pack sizes describe identical bit values, alternatively contiguous zeroes and contiguous ones, only the first value of the array needs to define the bit value of the first pack (see the red value in figure 3), and all subsequent bit values of the following packs can be determined by the decompression algorithm.
To give the possibility to data-processing tools to improve their memory footprint, we included special values in this encoding to mark the boundaries of AMR levels so that dataprocessing tools can perform on-the-fly decompression on a level-by-level basis for level-of-details approach. We used the 9 and 10 uint8 values as level boundary markers, two different values being required to mark level boundaries in case a bit flip occur (or not) at the end of an AMR level (see orange values in figure 4). In addition, we developed a special algorithm which is able to directly process the compressed stream to get the value at a given offset in the uncompressed stream without the need of uncompressing the entire AMR compressed array. We further reference this LightAMR grid boolean arrays lossless compression as base 52 continuous pack size encoding (CPS52). While this CPS52 encoding may seem peculiar compared to the majority of other RLE schemes from the literature, it yields very good performance. In the process of LightAMR data format standardisation though, this custom RLE would require to be updated to match more closely standard encoding schemes, prior to the LightAMR format standard release. Once the AMR grid description arrays are compressed, the next step in order to reduce the data volume is to compress the float (simple/double precision) scalar field data arrays.

LightAMR scalar field float data lossless compression
The development of our custom compression algorithm was motivated by the poor compression rates obtained with the deflate algorithm from the zlib 5 , about 3 to 5% with a compression speed about 50 MB/s. Even if higher compression rates could be achieved using the LZMA algorithm, it would require far more memory for a lower compression speed, such memoryhungry approach would not be viable in memory-bound massively parallel simulations as are many RAMSES runs. In addition, compression of float valued arrays using entropic encoding without prior knowledge of the data generally leads to poor performance due to its intrinsic homogeneous byte distribution unlike for example compressing a text where the distribution of byte (letters) is not homogeneous and therefore entropic compression is highly efficient.
In this section we present an adaptation of delta compression algorithm for lossless float data [17], [18] that we specifically designed for lightAMR octree structures. Delta compression is a commonly used compression method based on prediction functions. Indeed, this method is based on the choice of a mathematical predictor function and encode directly the prediction error. The compression process itself is performed by removing the leading zeroes in the encoding of the prediction error [17].
With this kind of approach, high compression ratios are reached when the predictor function is quite accurate, which is not easy to achieve when dealing with an AMR octree structure. Thus, we extended this delta compression algorithm to AMR grids by taking advantage of the fully threaded approach of the AMR strategy in RAMSES.
The only information our compression algorithm needs is the grid refinement array and a physical quantity field to compress (simple or double precision float). We propose to use an AMR-specific mathematical function to predict intensive variable scalar field values called the Parent-Child predictor (PCP) function. Indeed, the value of the parent cell is used as the predicted value for its child cells and so on from the top of the AMR tree down to the bottom, as represented on figure 5. The compression is then made octant by octant (if 3 dimensions are used). There are two mains steps in the data compression process: first, the delta compression itself with the leading zeroes removal and second, the encoding of the compressed scalar field on a bit-field .
As for the delta compression stage, the single (or double) precision parent float value of the physical field is mapped on unsigned integer of the same size operator to compute the differences between the mapped parent value and all the children ones, leading to a pack of prediction errors (s 1 , s 2 , s 3 , s 4 in the 2D example shown on figure 6). The next step is to compute an OR operation between all the s i values in order to determine how many leading zeroes could be removed on each s i value at most, which corresponds to the stage s 0 + s 1 + s 2 + s 3 in figure 6. In the example, we can remove a maximum of 5 leading zeroes on each s i value without losing information. At the encoding stage, we need to encode the residues of the prediction errors of all children cell values as well as the number of removed leading zeroes. On figure 7 we show the probability mass function (PMF) of the number of leading zeroes on cell basis and on pack basis on the double precision density scalar field of the Extreme-Horizon simulation (12800 density datasets). Since residue bits are close to being uniformly random, we cannot compress them well. If we compute the compression ratio in this example, without adding the encoding of the number of removed leading zeroes (nLZ), we obtain 1.303 on cell-by-cell basis and 1.237 on pack basis which is expected since more bits can be remove in the first scenario. Nevertheless, the number of removed leading zeroes for each cell or each pack needs to be encoded for the decompression stage. The shannon's entropy of the nLZ field on a cell-by-cell basis is 3.5 and 3.2 on a pack basis, this smaller value being also an expected result since we smooth the nLZ value in each pack. Nevertheless, the total number of nLZ values to encode in the cell-by-cell case is, in 3D, 8 times the total number of nLZ values to encode in the pack case. Therefore, even by adding another stage of compression of the nLZ values, we cannot hope to achieve a better compression ratio even if more bits can be removed from the data. If we take into account the encoding of both the nLZ values and the residues, the compression ratio drops to 1.217 for the the cell-by-cell case, while it only drops to 1.225 for the pack case. What is more, the cell-by-cell approach would require 8 times more computations of leading zeroes and an additional compression step to barely achieve a similar compression ratio. Thus, we chose to compute and encode the nLZ value once per pack of children as shown at the encoding stage on figure 6.
After computing the PMF of scalar field data, we notice that setting the number of removed leading zeroes on 4 bits by default is a good balance between compression speed and compression ratio. In the example of the Extreme-Horizon density field, increasing to 5 bits will only increase the compression ratio by 0.23 % but reduce the compression speed by almost 8%. Indeed, the nLZ computation is a key point of the algorithm performance and it requires more if-branches to compute a higher value. Finally, after the 4 bits used to encode the nLZ value, all residues s i are encoded, see figure 6. According to the number of bits used for the number of removed leading zeroes, the theoretical maximum compression ratio in dimension d is given by (E enc is the number of bits, 64 for double precision and 32 for single precision): Therefore, the maximum compression ratio achievable for single precision data is 1.829 and 1.293 for double precision with the default encoding. Nevertheless, we could expect this method to yield poor performance in the case of a difference of sign between a parent cell value and one of its child cells. Indeed, the very first bit of a negative value will be set to one thus no leading zeroes will be removed. This could be the case for example with a small velocity field component varying around zero. Finally, our delta compression algorithm works for intensive variables but could be extended to extensive variables if the predicted values for the child cells were taken as the value of the parent cell divided by the numbers of child cells. If the number of child cells was a multiple of 2, the mantissa of the prediction error would not be modified (only the exponent part would be bitwise right-shifted), which could still lead to a lossless compression algorithm. The RAMSES code using only intensive variable AMR scalar fields, this extension of the Parent-Child predictor function is beyond the scope of this work.
For the float data decompression stage, the AMR grid refinement array is used level-by-level to guide the delta decompression algorithm through the grid from top to bottom since the scalar fields values at level l are necessary as predictor values to decompress scalar field values at level l + 1.

Results and discussion
In this section we present the results of the lightAMR format conversion as well as tree pruning and compression results obtained on different real simulation snapshots : FRIG 6 , snapshot at 10.05 Myr and 1675 coarse steps, [19] composed of 4096 domains, ORION 7 snapshot at 1.26 Myr and 2125 coarse steps, [20] composed of 512 domains and Extreme-Horizon 8 , snapshot at Z ∼ 1 with 11866 coarse steps, [9], composed of 12800 domains.
Since, post-processing tools handle RAM-SES data on a per-domain basis and that lightAMR is a selfconsistent data model also on a per-domain basis, we consider a dataset as the data of one domain . For the tree pruning algorithm we show the minimum, maximum and average value for different metrics that are computed on a per-domain basis.

Tree pruning results
The table 2 gives the results of the tree pruning applied on the different datasets. The algorithm is 2.3 to 3.4 times more efficient on FRIG and ORION data than on Extreme-Horizon. FRIG and ORION runs are zoom simulations (nested grids with increasing spatial resolution) where most of the fine resolution grids can be found in the center or mid-plane of the simulation box, whereas Extreme-Horizon is a cosmological simulation where finest AMR levels are more homogeneously distributed in the box. These results are consistent with the fact that RAMSES is more memory efficient in cosmological runs than in zoom simulations. Nevertheless, the tree pruning gain is applied on the whole simulation data therefore, 11.42% on Extreme-Horizon means that the size is reduced by about 240 Gb. Moreover, tree pruning is done in a memory friendly way within a negligible time ( a few microseconds on the data of an MPI process ) . Therefore this step has a significant impact on data volume at almost no cost. Actually, the time spent doing tree pruning is correlated with the number of coarse cells and the position of those coarse cells in the array.
The efficiency of the tree pruning algorithm is far from intuitive in most cases (at least with RAMSES AMR grids). It depends on the maximum level of refinement, on the number of domains, on the geometrical convexity of the subdomains, on their surface to volume ratio, and is greatly impacted by the 2:1 balance constraint in RAMSES required by its internal (hydro) numerical schemes. The tree pruning rates can span a wide range within a factor 3.25, even for different domains of the same simulation snapshot (see

LightAMR grid boolean arrays compression results
We compared our compression speeds and ratios to commonly used and open source libraries : Zlib 9 , LZ4 10 and Snappy 11 . The compression ratio is simply the ratio between the uncompressed and the compressed buffer sizes. We run this compression benchmark on the Extreme-Horizon datasets and compare the compression speed and ratio of the cumulated grid refinement and ownership mask arrays after tree pruning, see table 3. Domain datasets are processed sequentially, the uncompressed and compressed sizes are accumulated as well as the times to compress them. The compression speeds and ratios are then computed and reported in table 3. While our single threaded algorithm is slower than the LZ4 package, we are able to achieve a higher compression ratio. The zlib at level 9 (best compression ratio) achieves a much higher compression ratio with the grid refinement array but at the cost of a significantly reduced speed. The compression speed of all those algorithms are linked to the layout of the uncompressed array, i.e. consecutive pack of zeroes or ones. With the RAMSES current parallelism (one MPI process per core) we barely have to compress arrays greater than a few millions of cells, even in the Extreme-Horizon use case. Therefore, a compression speed of 374 Mb/s is enough to compress the entire lightAMR arrays of the most populated MPI process of the Extreme-Horizon dataset in less than 8 ms (domain 427 with 3059513 cells after tree pruning). 9 Zlib website 10    The compression ratios and speeds are much higher with the ownership mask array because it contains very large packs of successive zeroes (cells that belong to the domain). This is due to the fact that boundaries generally happen on the border of a domain but also because the tree pruning algorithm removes ghost nodes. In addition, our algorithm performs much better at compressing such an array than the other tested libraries because of encoding large packs in base 52 produce a very small pattern.
One can notice that our encoding values are stored on 1 byte even if there are only ranging from 1 to 63, and therefore would require only 6 bits thus a potential gain of 25% in compressed size. In addition to that, computing the average Shannon entropy on 1 byte 'word' of all the AMR description array on an entire simulation like Extreme-Horizon leads to a value of 4. Which means that on average only 4 bits could be used to encode the information contained in one 'word'. The reason is simple, this is because all the number from 1 to 63 are not necessaryly used for the encoding. Actually, by computing exactly the number of bits that would be required to strictly encode bit to bit the information leads to a potential gain on average of 45% on Extreme-Horizon and thus to reduce the size on disk of the AMR data by 45% to less than 420 MB. We did not implemented such an optimization because it adds a step that will increase the time for compression for an AMR structure that is no longer significant compared to the size of the physical fields. And it will increase the complexity to rapidly retrieve AMR level seeds but also the direct access of cells states in the compressed data. Finally, in the case of Extreme-Horizon datasets, the AMR part of the data is now only representing 0.05% of the total simulation snapshot volume. Table 5 compare the com-pressed lightAMR format to the legacy RAMSES AMR format.
One important remark is that most data arrays stored in a RAMSES checkpoint are no longer stored in this new ligh-tAMR format, without any loss of information on the AMR grid structure itself. The large compression ratios presented in table 5 are mainly due to the switch to a boolean array description and the omission of this redundant meta-information, (e. g. neighbour cell indices, cell centers, parent cell indices) which can be easily and efficiently reconstructed on-the-fly by post-processing tools. The tree pruning step and the CPS52 encoding further increase this AMR file size compression ratio although in a less significant way.

Simulation
Ratio File size (MB) FRIG 1485.0 20 ORION 1500.0 6 ExH 583.3 840 Table 5: Chained ratio on RAMSES collections of datasets of tree pruning and CPS52 (left) and resulting lightAMR file size for the AMR grid description (right).

LightAMR scalar field float data compression results
Most of the data volume of an AMR simulation dataset is due to scalar fields represented as floating point data. The lossless compression algorithm takes the grid refinement array and an associated scalar field float data array as input, both after being processed by the tree pruning algorithm. We used a maximum number of removed leading zeroes of 15 (4 bits encoding) and the compression is single threaded. First, we compare in table 6 the compression performances of the different libraries on the whole density field (double precison, 12800 datasets) from the Extreme-Horizon datasets . We also integrate in the benchmark the Zfp library [21], which is more specific to floating point data compression. We show that our algorithm achieve a higher compression ratio on this dataset due to the use of the topology information from the AMR mesh. It was expected that dictionary compression approach used in LZ4, Snappy and Zlib cannot efficiently compress floating point data. Nevertheless, since the Zlib also use an entropy encoded, it is able to achieve a few percent of compression but with very slow speed compared to our algorithm.
On Extreme-Horizon, a single scalar quantity over one snapshot represents about 160 Gb (after tree pruning) thus a compression of 1.23 on the density field means a gain of about 30 Gb on disk. In table 7 we show the result of the PCP4 algorithm on the density, x-axis velocity, and thermal pressure for the different RAMSES datasets. We achieve a good compression ratio on the three tested scalar fields of about 1.19 in the worst case scenario and 1.23 on the best case. The compression speed is quite regular and above 1 GB/s, we noticed that this compression speed is linked with the CPU frequency thus with a higher CPU frequency we generally achieve a higher compression speed. For example, on an AMD Ryzen 7 3700X @  3.6 Ghz we were able to reach compression speed above 2.5 GB/s up to 3 GB/s on FRIG and ORION datasets.
We expected one current limitation of our algorithm for scalar fields with an average value of zero (e. g. velocity component). Indeed, if a parent value is positive and one of the children value is negative, then the very first bit of the integer representation will be set to 1 (sign bit) and therefore no delta compression could be achieved on that octant/quadrant. This limitation can be easily circumvented by shifting all values with a constant parameter. However, the results on table 7 show that the compression ratios for the velocity component are always at least as good as those of the density and pressure fields, highlighting the fact that the ocurrence of opposite sign between predicted value (parent cell value) and child value (child cell value) is rather uncommon in typical astrophysical use cases.
The main reason that explains our high compression speed is that we do not need to evaluate a mathematical function to predict each value but we merely use the parent value which is already available in memory. Moreover, this parent value is used as a predictor for all the child's cells (4 to 8 cells depending on the number of dimensions). Furthermore, based on the formula of the compression ratio, removing the bit of sign plus the entire exponent on the 64 bits integer representation of a double precision float for a pack of cells leads to a theoretical ratio of 1.219. We generally achieve a ratio higher than this threshold according to the results in table 7, which means that almost every child cell value lies within a factor two of the parent cell value, which demonstrates the interest of using our predictor function. Since the basic idea of AMR is to refine a cell to avoid high variation of the physical field between adjacent cells, it was expected that in most cases we achieve at least this ratio. Nevertheless in the density field of FRIG, this is not the case mainly because a lot of shocks (turbulent interstellar medium modeling) are present in the simulation, a large number of steep density discontinuities can lead to lower compression ratio.
Since the float compression data algorithm process the AMR tree from the top to the bottom, from coarse cells to leaves, we can not access directly the value of a cell at a level l without first decompressing the tree from level 0 up to level l. There-   Since in RAMSES one can downcast some physical quantities to single precision in order to reduce data volume by a factor 2, we give in table 9 the compression speeds and ratios for single precision data. If the entire exponent part (8 bits) of a single precision float (mapped on a 32 bits unsigned integer) is removed during the delta compression step, the theoretical compression ratio is 1.362, while it was only 1.219 for a double precision float. As shown previously with double precision float compression, we are still able to remove in most of the datasets the entire sign plus exponent part of the bit representation in the single precision float compression. Since the sign and exponent represent a bigger part of the 32 bits representation, it naturally leads to higher compression ratios. Therefore, downcasting double precision physical scalar fields within RAMSES can not only yield a reduction factor of 2 (due to down-casting) but also an even more efficient compression using our PCP algorithm. Consequently, down-casting from double to single precision can be of great interest in order to significantly reduce data volumes if double precision is not mandatory for the post-processing workflow.  The PCP4 algorithm is based on parent value as predictor value, this algorithm cannot be used if coarse node value are not available. Nevertheless, the basic idea of the lightAMR data model is to be a post-processing dedicated data model. AMR postprocessing tools heavily use LOD approach which requires node value in order to be efficient. Therefore, if no node value is provided, then the whole dataset must be loaded to compute node values through up-sampling, such an operation will introduce very large overhead even for simple task such as low resolution images of the simulation box. The typical compression ratios that we obtain more than compensate the 12.5 % overhead required to store node values (in 3D cell-based AMR octree).

Conclusion
Adaptive Mesh Refinement suffers from a lack of standardized format to describe the grid structure and therefore AMR codes generally use their own format or eventually the unstructured grid to describe the mesh. Unfortunately, using the unstructured grid description is highly inefficient with AMR meshes and produce very large file even for relatively small meshes. In this paper, we presented an extremely compact format designed for AMR meshes, called lightAMR, for which we developed a set of dedicated data reduction and compression algorithms. To minimize data redundancy in RAMSES AMR grids, we developed a tree pruning algorithm able to remove between 11% and 38 % of cells in various astrophysical datasets. The large discrepancy of the obtained values is mainly due to the diversity of adaptive refinement layouts in our selected dataset simulation boxes. On top of that, we implemented a lossless and memory friendly compression algorithm specific to the lightAMR format that compresses the grid refinement/ownership mask arrays using base 52 continuous pack size, called CPS52. By switching to a boolean array description of the AMR grid in the lightAMR format, by omitting unnecessary meta-information and by chaining the tree pruning and the CPS52 algorithms on the tested RAMSES AMR datasets, we obtain AMR mesh file size compression ratios ranging from 583.3 to 1500.0, thus reducing the AMR grid description to a negligible part of the overall data volume.
In addition, a lossless, memory friendly and efficient floating point data compression algorithm, called PCP, designed for lightAMR physical scalar fields achieve a high compression ratio and speed compared to commonly used and open source libraries by taking advantage of the topology of the AMR structure. For double precision, we achieve a compression speed ranging from 960 MB/s to 1275 MB/s and a ratio between 1.17 to 1.26, on RAMSES datasets, with a sequential version of the algorithm and a Intel Xeon Gold 5118 @ 2.30 GHz. While down-casting physical scalar fields to single precision naturally reduces the size by a factor 2, it also increases the compression ratios that range from 1.30 to 1.54 with similar compression speeds.
Finally, for the tested astrophysical datasets, the use of the lightAMR format leads to overall data reduction percentages of 62.26 % on FRIG, 49.64 % on ORION and 37% on Extreme-Horizon with no loss of information. If a RAMSES user is ready to sacrifice the precision of the scalar quantity fields and lets the lightAMR format downcast float data to single precision, then he can expect even more significant data volume reduction percentages (83.2% on FRIG, 81.3 % on ORION, 75.8 % on Extreme-Horizon). In the context of the three astrophysical simulation projects tested in this work, all the while using the same disk storage capacity, far more ambitious data output policies could have been considered, with data output done at higher frequencies (e.g. x5 on ORION, x6 on FRIG, x4 on Extreme-Horizon). The lightAMR format, as well as compression algorithms, are compatible by design with level-of-details (LOD) approach which is a very important technique used by post-processing tools to improve performance for both visualization and data analysis.  Table 10: Recap of data volume reduction and compression on RAMSES simulation datasets. LightAMR data takes into account the tree pruning algorithm. The first two lines gives the total number of cells before and after tree pruning. The AMR grid data volume using both the tree pruning and the CPS52 compression algorithms and the HYDRO data volume is reduced using the tree pruning and the PCP4 compression algorithms. The last two lines show the data volume reduction when float data are down-casted to single precision (SP).
Nevertheless, we must insist on the fact that the lighAMR data structure is designed to significantly reduce the data stor-age impact of RAMSES (or any other fully threaded octree code) simulations as presented in the paper and must be seen as a "storage data structure". For analysis purposes, a "computationalfriendly" data strucuture, [23], should be constructed from the lightAMR data in order to achieve good performance. In our analysis workflow, for example, we use vectors or hashmaps of leave cells only, allowing us to achieve high performance data analysis with an hybrid multiprocessing / multithreading parallelism but this discussion is beyond the scope of this article.

Future work
In the context of the standardisation process of the data format, the custom CPS52 run-length encoder will be updated in the near future to match more closely other mainstream runlength encoding schemes found in the literature, prior to the first release of LightAMR format. While the detailed PCP4 compression is lossless, we will explore in future work lossy compression as a compromise to down-casting physical quantities to single precision in order to further reduce I/O volume. In order to further reduce the data volume for post-processing and since the lightAMR data model is self-consistent and selfdescribing even over one domain, we will explore the concatenation of multiple ligthAMR data from adjacent domains. Finally, we plan to use this updated version of RAMSES with the integrated lightAMR format and parallel I/O with Hercule library for a large and ambitious astrophysical simulation run on a pre-exascale architecture.