L Rezeau

G Belmont

R Manuzzo

N Aunai

J Dargent

Analyzing the magnetopause internal structure: new possibilities offered by MMS tested in a case study

The internal structure of the magnetopause is investigated, using new analysis tools allowed by the high performance MMS instruments.

• In a case study, the observed boundary is shown to be non planar and non stationary which makes it necessary to perform a local study of the internal structure of the boundary. • Thanks to this local analysis, quasi 1D thin sub layers are identified separated by regions that are mainly 2D.

Introduction

The magnetopause boundary separates two magnetized plasmas of different origins: the solar wind and the magnetosphere. Its existence is due to the frozen-in property that prevails at large scale and which would fully prevent the two plasmas to reconnect if it was valid always and everywhere. As the magnetopause is accessible to in-situ spacecraft measurements, it provides a unique occasion to study the internal structure of such a boundary and understand how the two plasmas interpenetrate each other via the kinetic effects. However, this study is made difficult by the fact that the boundary is always perturbed by non stationary effects, due to the non stationary incident solar wind and/or to surface wave instabilities such as tearing and Kelvin-Helmholtz instabilities [START_REF] Chen | Tearing instability, Kelvin-Helmholtz instability, and magnetic reconnection[END_REF], [START_REF] Faganello | Competing mechanisms of plasma transport in inhomogeneous configurations with velocity shear: The solarwind interaction with earth's magnetosphere[END_REF]). It is worth noticing that, if purely planar and stationary, the magnetopause layer should obey the classical theory of discontinuities [START_REF] Belmont | Collisionless plasmas in astrophysics[END_REF], i.e. be purely tangential (B N = 0) or, if not, either purely rotational or purely compressional. This is in contradiction with observations since compressional and rotational variations are always observed in a close vicinity of each other in the magnetopause layer, often mixed but with sometimes a clear separation between both [Dorville et al., 2014a]. Thanks to its unprecedented high quality and high time resolution experiments, the MMS spacecraft [START_REF] Pollock | Fast plasma investigation for magnetospheric multiscale[END_REF], [START_REF] Russell | The magnetospheric multiscale magnetometers[END_REF]) nowadays allow significant advances in the study of the internal structure of the magnetopause layer. This paper shows the new methods that can be used for that purpose. October 16th, 2015 was a day with multiple magnetopause crossings by MMS. Fig. 1 shows that it is due to the fact that the orbit of the spacecraft grazes the magnetopause for about 4 hours between 09:00 and 13:00 UTC. The expected position of the magnetopause is calculated with the Shue et al model ( [START_REF] Shue | A new functional form to study the solar wind control of the magnetopause size and shape[END_REF]) using ACE data [START_REF] Stone | The Advanced Composition Explorer[END_REF]. The figure evidences that many crossings are expected to happen. This is what is observed and these multiple crossings can be expected to be complex, with possible back and forth motions and partial penetration in the current layer. We choose to study the crossing around 13:06 (which is shown by a red arrow) because this period has already been studied by [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF], [START_REF] Torbert | Estimates of terms in ohm's law during an encounter with an electron diffusion region[END_REF], and Le [START_REF] Contel | Whistler mode waves and hall fields detected by mms during a dayside magnetopause crossing[END_REF], with a main emphasis put on its relationship with the reconnection event identified at 13:07. Fig. 2 displays the magnetic field measured by the MMS magnetometers [START_REF] Russell | The magnetospheric multiscale magnetometers[END_REF] during a 1 minute interval around the crossing investigated. In this figure as in all the others unless specified, the times are counted for convenience from t 0 =13:05:30. The magnetic field is smoothed using a gaussian filter, with a standard deviation of the gaussian kernel equal to 70 points, which makes an effective smoothing window of about 1.6s. All the data used in the study are resampled to the magnetic field sampling time and then smoothed in the same way as the magnetic field.

One can see that the crossing is complex. The spacecraft come from a clear magnetospheric field at the beginning of the interval (B z ≈ 35nT); a reversal is seen around t ≈ 15s, showing the crossing of the main magnetopause current layer; the magnetic field is not completely stationary afterwards, which can be interpreted, as done by Torbert et al, by the fact that the spacecraft do not progress further in the normal direction with respect to the magnetopause, so remaining inside it ("stagnation"), with even a backward motion around t = 28s. Fig. 3 summarizes the evolution of the main physical parameters during the interval under study, where it can be seen that the region where the plasma properties change is not identical to the magnetic field reversal region but is close to the first part of it, and slightly before.

The magnetopause is non-stationary and non-planar

Comparison of normals

The most common method to analyze a magnetopause crossing is the Minimum Variance Analysis (MVA), which has been introduced with the first measurements of the magnetic field in space ( [START_REF] Sonnerup | Magnetopause structure and attitude from explorer 12 observations[END_REF], [START_REF] Sonnerup | Minimum and Maximum Variance Analysis[END_REF]). It is based on the assumption that the boundary is perfectly 1D, i.e. that all isosurfaces are parallel planes, and it provides a single boundary normal based on the magnetic field measurements across the "whole crossing". Years of study of experimental results have shown that this assumption is acceptable as long as sufficiently large scales are considered and ultimately amount to finding out the normal of the magnetopause boundary itself and compare it to a model, e.g. [START_REF] Shue | A new functional form to study the solar wind control of the magnetopause size and shape[END_REF]. But they have also shown that the magnetopause itself has an internal structure which can be complex ([Dorville et al., 2014a], [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF]).

MVA relies on the Maxwell equation ∇ • B = 0, and on the constancy of the normal component that follows from it for a strictly 1D geometry. This property is sufficient to determine the normal direction as long as this component is the only that does not vary, i.e. when the B T tangential hodogram has a certain curvature: otherwise, two components are constant and B N =cst is not a sufficient condition to identify the normal direction (this excludes the coplanar case of shocks). If the magnetopause conformed to the simple classical image of a boundary made of a monotonous ramp connecting two homogeneous regions, the strict B N conservation would be valid on any interval, whatever the number of points. The existence of different sub layers that can move with respect to each other would not invalidate this property, at the condition that these sub layers are all planar and strictly parallel to each other. The existence of non stationarity should not bring difficul-ties either, at the condition that the boundary remains strictly planar everywhere and that its normal direction does not vary in time. The main difficulties therefore come from the departures from planarity and from the absence of time stationarity of the normal direction. Such departures are likely to occur often at the magnetopause, even if only due to the small scales waves and turbulence that are always present. To fix this difficulty, MVA is usually used on a statistical basis and applied over a sufficiently long interval between two points around the crossing, one considered as assuredly in the magnetosphere and one as assuredly in the magnetosheath. This actually transforms the condition that B N is constant into the condition that its variance is less than the variance of the other components. A necessary condition for applying safely this condition is that the ratio between the minimum and intermediate variances is sufficiently small. Another condition that should be checked is that these two variances are really characteristic of the large scale variation related with the current layer under study and not mainly due to the parasitic small scale turbulence. When these conditions are not fulfilled, the result actually depends on the position and the size of the "global" interval chosen. The stability of the result is sometimes tested a posteriori, by checking the variations of the observed B N and by using nested intervals (see for instance [START_REF] Zhang | Neutral sheet normal direction determination[END_REF]). When contradictions occur in one of these two tests, the results are rejected, under the assumption that the real local normal should not depend on time inside the crossing. Beyond this constraint of a strictly constant normal direction, MVA also suffers from another limitation that prevents people from using it on short intervals and therefore analyzing the sub-structure of the layer: the interval used must be long enough to evidence the curvature of the B T tangential hodogram. Any variation obviously tends toward a straight line when the interval duration decreases, so increasing the inaccuracy of the result in the M-N plane.

These limitations encouraged scientists to develop more elaborate methods (a review can be found in [START_REF] Haaland | Four-spacecraft determination of magnetopause orientation, motion and thickness: comparison with results from single-spacecraft methods[END_REF]). They are not all used nowadays, probably because they require more high time resolution data and are more difficult to apply than MVA. Let us cite in particular the different GRA methods (Generic Residue Analysis [START_REF] Sonnerup | Orientation and motion of a plasma discontinuity from single-spacecraft measurements: Generic residue analysis of Cluster data[END_REF]). These are generalizations of MVA to other parameters than just B. Although generally more efficient than MVA, these methods rely on conservation laws (fields and plasma) that require also planarity (1D variations) to be valid. They therefore suffer from most of the limitations of MVA for investigating sub-layers. In addition, they require stationarity (∂ t = 0). The BV method [Dorville et al., 2014b] mixes magnetic field (B) and velocity (V) data and is based on different grounds but still in the same "global layer" spirit. It has been shown to give accurate normal determinations in a statistical study [START_REF] Dorville | Magnetopause orientation: Comparison between generic residue analysis and bv method[END_REF]. Nevertheless it is not either perfectly suited for analyzing intervals much shorter than the global width of the current layer (in spite of the excellent time resolution of the MMS data). In any case, all the methods mentioned here assume the boundary is locally a plane. This assumption may be questionable due to local deformations of the surface, such as surface waves. Confirmation is given by all the numerical simulations of reconnection or Kelvin-Helmholtz instability show [START_REF] Aunai | Orientation of the x-line in asymmetric magnetic reconnection[END_REF], [START_REF] Chen | Tearing instability, Kelvin-Helmholtz instability, and magnetic reconnection[END_REF], [START_REF] Dargent | Kinetic simulation of asymmetric magnetic reconnection with cold ions[END_REF], [START_REF] Miura | Nonlocal stability analysis of the mhd kelvinhelmholtz instability in a compressible plasma[END_REF]) and as some experimental observations [START_REF] Blagau | A new technique for determining orientation and motion of a 2-D, non-planar magnetopause[END_REF].

For the crossing investigated in the present paper, MVA has been first applied on the global interval. It shows that the three eigenvalues are not well separated, the maximum variance being clearly larger than the two others, but these two others being rather similar (ratio 1.9). This means that the normal might not be precisely determined. Nevertheless, we obtain N MV A = [0.811, 0.536, -0.234] which is close to the normal obtained with the [START_REF] Shue | A new functional form to study the solar wind control of the magnetopause size and shape[END_REF] model which is N Shue = [0.854, 0.519, -0.043]. The angle between the two normals is 11°indicating that in this case the "global" magnetopause is probably not far from the standard paraboloid shape assumed by Shue et al. As MVA, as we use it, is a single-spacecraft technique,one can compare the MVA normals derived from the data on each of the four spacecraft. As they are actually very near, they measure very similar magnetic fields and the angle between each normal and the average normal is indeed less than 1°.

Looking at the magnetic data in Fig. 2, the global crossing can be guessed to consist of a first current layer between, typically, t = 10s and 20s, followed by a backward motion later, with only a partial entrance in the magnetopause between t = 25s and 30s. For confirming or disproving such a guess, one has to investigate the internal structure of the magnetopause layer in more details and look for possible sub-structures. For this purpose, let us first compute MVA on shorter intervals. Between t = 10s and 20s, we obtain (on MMS1) N MV A = [0.591, -0.591, -0.548], which is very different from the previous normal, the angle between both being 73°. Let us note that changing slightly the choice of the beginning and ending times of this MVA interval does not change much this conclusion. As the ratio between minimum and intermediate eigenvalues is again not much larger then 1 (2.6), MVA is quite questionable and one can wonder whether this determination is just erroneous or if such a large difference can actually exist between the local and the global normals. Taking advantage that, beyond B, all the other physical parameters are measured at the same time, it is possible to use the particle data [START_REF] Pollock | Fast plasma investigation for magnetospheric multiscale[END_REF] to analyze the crossing with the BV technique [Dorville et al., 2014b]. The hodogram (Fig. 4) is almost a straight line, without a clear curvature, but this does not prevent the BV method from working, the fit of this hodogram by a very elongated ellipse remaining quite acceptable. The BV program automatically determines the optimum interval for its fitting procedure, which is between, unsurprisingly, t =14s and 18s. The normal obtained is then (on MMS1): N BV = [0.838, 0.506, -0.205], which is only 9°from the Shue et al normal. This result is much more likely than the MVA one.

Thickness of the magnetopause

A possible byproduct of the BV method is an estimation of the thickness of the current layer of the magnetopause and of its normal velocity, but it is worth noticing that these estimations have to be taken with caution. The BV program provides, in its present version an estimated thickness of 30 km on MMS1 and MMS2 and 40 km on MMS3 and MMS4, which is smaller than the thermal ion Larmor radii (which vary from ≈ 140 km in the magnetosphere to ≈ 110 km in the magnetosheath). It also provides an estimated normal velocity of 8 km.s -1 for MMS1 and MMS2 and 10 km.s -1 for MMS 3 and MMS4, which is much smaller than the normal Alfvén velocity (36 to 170 km.s -1 ). These results being noticeably smaller than the values commonly observed, we have used other methods to check them. These methods provide more likely results of about 200 km for the thickness and 50 km.s -1 for the normal velocity.

The first calculation is the same as done in the BV method, but also similar to those used in [START_REF] Paschmann | The magnetopause and boundary layer for small magnetic shear -Convection electric fields and reconnection[END_REF] and [START_REF] De Keyser | Trying to bring the magnetopause to a standstill[END_REF], which consists in integrating the normal ion velocity V in over time to obtain the abscissa x(t), but using a different normal which is likely to be more precise (see in further sections how we have obtained this normal). The second calculation makes use of the four-spacecraft gradient determination. The abscissa along the normal is obtained by integrating the quantity δx = Y /Y , where Y is a scalar variable and where Y represents the projection of ∇Y on the normal direction (the normal direction being determined in the same way as above). The spatial derivatives in the different directions are estimated by linear interpolations from the multi-point measurements (here 4 spacecraft). This can be done by methods similar to the well-known "curlometer", which is very often used to calculate the electric current density [START_REF] Chanteur | Spatial Interpolation for Four Spacecraft: Theory[END_REF]. We have taken here Y = B L , which is the component of B that varies most during the crossing. Fig. 5 shows the comparison between the two results. Both results look quite compatible during the crossing of the main current layer and lead to the same value of ≈ 200 km for its thickness. This similarity validates the hypothesis which is done in the BV method that the flow through the boundary is negligible. Nevertheless, the two results clearly depart at later times. This is due to a very strong dependence of the result, with the BV method, on the quality of the normal determination [Dorville et al., 2014b]. A small uncertainty in the normal direction determination can draw a large variation of the V n component because the tangential component of the velocity is much larger than the normal one (see Fig. 3). With a magnitude of the velocity of about ≈ 300 km.s -1 , an uncertainty of 10°on the normal direction corresponds to an uncertainty of ≈ 50 km.s -1 for the normal velocity, and an uncertainty of about ≈ 200 km for the thickness. It is so quite understandable that, with a normal valid in the 14-18s interval, the inaccuracy increases very fast at later times where this normal is no more valid. The method based on gradients does not present this difficulty: it is much less sensitive to the accuracy of the normal determination. Nevertheless, we had also to add a caution to make it work correctly: because of various small accuracy issues, the denominator Y may cancel at a time slightly different from the numerator, which results in short divergences in the result and jumps in the x(t) curve. This has to be corrected by adding adequate small shifts in the denominator. In addition, Fig. 5 clearly gives the confirmation that the spacecraft is going backward inside the magnetopause around t = 27s, as was guessed before. Due to its importance, this technique is under review for further improvements and will be applied to other cases in next studies. In Fig. 6, we have plotted the projection of the ion velocity along the normal obtained by BV, together with the density profile. This evidences an internal structure inside the magnetopause. Two main parts can be observed in the interval t = 14 -18s, where the main plasma gradients are located and which is emphasized by a thick line: in interval (a) a sharp density gradient, with an almost constant V n , followed in interval (b) by a smoother gradient with a normal velocity close to zero. This is in agreement with the sketch drawn in Fig. 3 of [START_REF] Burch | Electron-scale measurements of magnetic reconnection in space[END_REF] which is a possible interpretation of this crossing (although assuming a stationary boundary): a rather straight crossing, followed by a stagnation of the spacecraft inside the boundary. This is confirmed by the observation of energetic ions continuously after 13:05:42 [START_REF] Contel | Whistler mode waves and hall fields detected by mms during a dayside magnetopause crossing[END_REF]. Out of the central interval t = 14 -18s, the curve V n (t) is plotted with a dashed line, to warn the reader that the projection of the velocity is obtained using the BV normal based on this interval and that the validity of this projection, even if correct in the magnetic ramp itself, remains questionable outside of it.

Non-stationarity

Using timing methods is another classical way for getting information on the boundary properties from multi-spacecraft measurements. We tested CVA (Constant Velocity Analysis), which assumes the boundary is a planar structure encountered by the 4 spacecraft with a constant velocity [Sonnerup et al., 2008a], [Sonnerup et al., 2008b]. As in any other timing method, the analysis is based on the knowledge of the positions of the spacecraft and the measurements of the delays between the signatures of the crossing seen by the four spacecraft. As shown in Fig. 7, these delays are very short with respect to the parasitic variations due to the intrinsic non stationarities, in particular waves and turbulence. If the boundary was stationary, we should find a constant delay between the fields observed by MMS1 and MMS4. On the contrary, it is obvious that the dispersion of the points is not negligible at all with respect to the delay itself. It is worth noticing that we have plotted here the B z component, which is the component that varies most, and for the MMS1-MMS4 pair, for which the delay is maximum. The situation is worst when using the other components and the other spacecraft pairs. This results in a very inaccurate determination of the delays and therefore in a bad determination of the normal direction. The first conclusion is therefore that, in this case, the CVA method cannot be used without much caution.

Looking at Fig. 7, we can also derive some hints on the non stationarity of the boundary at different scales. In the beginning of the crossing there are oscillations, evoking the presence of waves, superimposed to the magnetopause variation. This induces variations of the delay on the top of the figure. But there is also a large-scale variation of the delay: on the top of the figure (beginning of the crossing) its mean value is about -0.07s and afterward it goes to -0.15s: the delay is not constant through the crossing. Similar conclusions are obtained with the two other spacecraft. Using an averaging of the delays, one could interpret the large-scale variation as a constant acceleration of the boundary, which would help improving this result [START_REF] Dunlop | Four-point cluster application of magnetic field analysis tools: The discontinuity analyzer[END_REF]. Results of other timing methods, such as CTA (Constant Thickness Analysis) are not presented here, but the same difficulty (small delays with respect to the intrinsic fluctuations) would lead, on this example, to the same difficulties. time (s) from 13:05:30 delays (s) between the magnetic fields The conclusion of these observations is that the magnetic field is not stationary during the crossing by the four MMS spacecraft and therefore the boundary is not the planar stationary discontinuity which is the most simple model for the magnetopause. It is necessary to investigate in more details the geometry and behaviour of the magnetopause.

Internal structure: departures from planarity

When analyzing a boundary crossing, one most often assumes that this boundary is 1D, i.e. that all parameters vary only in one direction, which is its normal. When this hypothesis of planarity is fully verified, the normal component B n of the magnetic field is strictly constant and this property is used in MVA method to determine a single "global" normal direction (if no other B component is constant in the interval). Nevertheless, when the boundary is shaken by some non stationary effect (either due to varying incident conditions or due to surface instability such as tearing mode or Kelvin Helmholtz), it generally does not remain fully 1D. Such departures to planarity can easily be observed in numerical simulations of reconnection (see for instance [START_REF] Dargent | Kinetic simulation of asymmetric magnetic reconnection with cold ions[END_REF], which will be used afterwards in the paper) or, less easily, it can be guessed from data (see the magnetopause reconstructions in [START_REF] Hasegawa | Optimal reconstruction of magnetopause structures from Cluster data[END_REF], De [START_REF] De Keyser | Empirical Reconstruction[END_REF]). These departures result in the fact that MVA is not suitable to this case and the meaning a global normal direction becomes unclear. One way for dealing with these cases is to try to determine, when possible, a "local normal", possibly varying along the crossing, instead of a single "global" one.

Local Normal Analysis

We introduce here a new method, that we call LNA (Local Normal Analysis), based on the independent measurements of B (from field data) and j (from particle data), and which allows determining a normal that can vary along the crossing. Mathematically speaking, a local normal direction can be defined wherever all plasma parameters depend on space only through a single scalar function s(x, y, z) of the three coordinates. This ensures that the gradients of all parameters are parallel to each other at any point, this common "normal" direction possibly depending on the point considered. The direction N is given by:

N = ∇s |∇s| (1)
In a cylindrical geometry for instance, all quantities depend on space only through the radius r, so that all gradients are everywhere parallel to the radial direction. Of course, this direction is variable from one point to another in the azimuthal direction.

For any vectorial field U verifying this property, one can write the curl as

∇ × U = ∇s × d s U =| ∇s | N × d s U (2)
where d s U stands for the derivative of U with respect to s. Therefore when it is applied to the magnetic field it shows that the current density is perpendicular to the normal (neglecting the displacement current). When applied to the electric field, it shows that ∂ t B is perpendicular to the normal, using Maxwell-Faraday equation. A simple cross product between these two vectors is then a priori sufficient to provide the normal direction

N = j × ∂ t B |j × ∂ t B| (3)
When both parameters j and B are independently determined with a sufficient accuracy, this expression can provide a simple and efficient way for determining the local normal N at each time and for a single spacecraft. It is worth noticing that this method does not rely on ∇ • B = 0 and thus on the fact that one component (and only one) is constant: it is therefore not limited to sufficiently rotational cases. For the first time in space history, MMS provides independent -and generally reliable-measurements for j and B [START_REF] Torbert | Estimates of terms in ohm's law during an encounter with an electron diffusion region[END_REF], since we can compute a high resolution current density from the particle data [START_REF] Pollock | Fast plasma investigation for magnetospheric multiscale[END_REF]. On previous space observations we used to work only with current density obtained from the magnetic field, with the well-known curlometer technique, because the particle instruments had neither the necessary accuracy nor the necessary time resolution to do it. On MMS it has been shown that both calculations of the current show a global fairly good agreement (see [START_REF] Contel | Whistler mode waves and hall fields detected by mms during a dayside magnetopause crossing[END_REF] who computed the currents for the same time period).

It is worth noticing that this new method has to be scale dependent: in the present program, this dependence is crudely controlled by the way the variables are smoothed before use. Since the method relies on time derivatives, this smoothing has an important role in the result. Here, the components of the magnetic field are smoothed with a local cubic fit, which is convenient for getting the time derivatives analytically (The smoothing is performed on the same timescale as the previous gaussian filtering). Going to large-scale smoothing should allow retrieving the classical notion of global normal. On the contrary, going to very short scale smoothing would provide the wave vectors of the different waves encountered (which can be considered as "parasitic" for the present kind of study). This step could be improved in the future (by using for instance a Fourier filtering instead of a smoothing). Fig. 8 shows what the results look like when running the "Local Normal Analysis" (LNA) method on the case presented in Fig. 2 without further precaution. The data have been smoothed over 1.6 seconds (the global interval being of 1min). This time scale is a good compromise for this case: it is significantly shorter than the global crossing time (so giving access to the internal structure), and long enough to get rid of most high frequency turbulence. One can see that this figure appears almost unintelligible in these conditions: apart from a short period about t = 15s where the normal appears relatively stable (and where its direction will be confirmed by another method hereafter), it appears highly fluctuating and apparently random. The reason can easily be understood: the method provides the local normal under the hypothesis that this normal exists, i.e. that the variations are locally 1D. As, at this stage, there is no test of this hypothesis, one gets a result everywhere, even where it is not verified and where the result is thus meaningless. An additional test of locally 1D variations is therefore necessary to make the LNA method complete. It will be the subject of the next sections.

Test of the local planarity

The best test for determining the dimensionality of observed variations demands multi-point measurements. It has been proposed by [START_REF] Shi | Dimensional analysis of observed structures using multipoint magnetic field measurements: Application to Cluster[END_REF] for Cluster data. This method, called MDD (Minimum Directional Derivative) analysis makes use of magnetic field data, although it is not based on specific properties of this field. It actually has been little used with Cluster, most of the authors preferring to stay in the purely 1D hypothesis and the simple notion of a global normal supposed to be determined by MVA. But it is nowadays attracting increasing interest for analyzing the MMS data (see for instance [START_REF] Chen | Electron diffusion region during magnetopause reconnection with an intermediate guide field: Magnetospheric multiscale observations[END_REF]) because of the short separation between spacecraft that allows a better determination of the local gradients. In a recent paper, [START_REF] Denton | Motion of the mms spacecraft relative to the magnetic reconnection structure observed on 16 october 2015 at 1307 ut[END_REF] have even applied this MDD method on a magnetopause crossing in the same global interval shown in Fig. 1 as the crossing analyzed here, but a bit later.

The MDD method consists in diagonalizing the matrix L = G • G T , where G = ∇B and the superscript T indicates matrix transposition and where the spatial derivatives are computed as explained before. The largest eigenvalue λ 1 corresponds to the largest deriva-tive for the ensemble of the B components. When this eigenvalue is much larger than the two other eigenvalues, it means that all B components vary in one single direction, which is given by the corresponding eigenvector α 1 , i.e. that it is 1D, with the normal direction N = α 1 . When the two largest values λ 1 and λ 2 have the same order of magnitude, while the third one λ 3 is much smaller, it means that the problem is 2D, the variations occurring in the plane (α 1 , α 2 ), α 3 so being the direction of invariance. When the three eigenvalues have the same order of magnitude, it means that the B variations are fully 3D. A modified MDD method has been proposed by [START_REF] Denton | Test of methods to infer the magnetic reconnection geometry from spacecraft data[END_REF] (see also a test in simulation in [START_REF] Denton | Test of Shi et al. method to infer the magnetic reconnection geometry from spacecraft data: MHD simulation with guide field and antiparallel kinetic simulation[END_REF]) to avoid the effects of possible offsets and calibration errors in the data. These errors might have a noticeable impact when the method is used to compute the velocity of a structure [START_REF] Denton | Test of methods to infer the magnetic reconnection geometry from spacecraft data[END_REF]) but, as it is not what we do here, we use only the original version of MDD in the present paper. Nevertheless, this point of view may have to be reconsidered for the generalized MDD method that we propose hereafter because such errors have certainly a much larger effect when using the electric field data than with the only magnetic field ones.

In order to visualize more easily the effective dimensionality of the variations, we have introduced three parameters, which can be used as proxies:

D 1 = λ 1 -λ 2 λ 1 (4) D 2 = λ 2 -λ 3 λ 1 (5) D 3 = λ 3 λ 1 (6)
These three parameters vary between 0 and 1 and their sum is equal to 1. For D 1 = 1 and D 2 = D 3 = 0, variation happens only in one direction: the geometry can be told "purely 1D variation". For D 2 = 1 and D 1 = D 3 = 0, the amplitudes of the variations are equal in two directions: it is what we call the case "purely 2D". For D 3 = 1 and D 1 = D 2 = 0, the amplitudes of the variations are equal in the three directions: it is what we call it "purely 3D". Of course, all intermediate situations are possible. Let us consider, for instance, a flux rope with λ 1 = 5, λ 2 = 1 and λ 1 = 0.1, which gives the dimensions D 1 = 0.8, D 2 = 0.18 and D 3 = 0.02. The structure has a slightly 2D character since D 2 is not negligible, but D 1 > D 2 indicates that the tube is strongly flattened in one direction: this makes the transition between 2D (circular tube) and 1D (tube infinitely flattened). Such structures have been observed and studied by [START_REF] Shi | Spatial structures of magnetic depression in the earth's high-altitude cusp: Cluster multipoint ob-servations[END_REF] and [START_REF] Shi | Solar wind entry into the high-latitude terrestrial magnetosphere during geomagnetically quiet times[END_REF] on Cluster and [START_REF] Yao | Observations of kinetic-size magnetic holes in the magnetosheath[END_REF] on MMS.

When applying the MDD Analysis to the interval under study, the three eigenvalues obtained are quite similar to those of the Fig. 1 of [START_REF] Denton | Motion of the mms spacecraft relative to the magnetic reconnection structure observed on 16 october 2015 at 1307 ut[END_REF]. These results are plotted in Fig. 9 using the three D i parameters. It must be kept in mind that the D i coefficients deriving from MDD give a local measurement of the dimensionality at the scale which has been selected by the smoothing. Our data have been smoothed on 1.6 s, therefore the wave structures superimposed on the magnetopause crossing are mostly removed. It can be observed that the 1D variations are generally dominant but that 2D and 3D variations are also present in the interval. It is worth noticing that, in the regions of 2D variations, the direction of invariance α 3 is determined by the MDD method, which may be an important information for numerical modeling purposes.

In the regions where D 1 ≈ 1, the normal can be determined by N M DD = α 1 . In Fig. 10, the angular distance of this MDD normal with the reference N Shue normal is plotted, for the regions where D 1 > 0.9 (thin line) and for D 1 > 0.98 (thick line). An additional caution has been taken in this figure: we have discarded the regions where there are no significant magnetic field variations (|∂ t (B)| 2 less than 1/10 of its maximum value) because we are not interested in the direction of the gradients for these small variations: they are more likely related to wave and turbulence rather than to the large-scale current layers. In the remaining regions, the results of our LNA have been over-plotted for comparison (in blue). One can observe that, as expected, the results obtained by the two methods are generally close to each other when D 1 ≈ 1, and that they diverge from each other for smaller values of D 1 . For the sake of clarity, we have isolated the two intervals, limited by dashed lines in the figure, where D 1 > 0.98 and which are long enough: interval 1 from 13.8 to 16.8, and interval 2 from 27.4 to 28.4. If we compute the averaged normals on these intervals, we find that the two normals make a 4°angle in interval 1 and 7°angle in interval 2. Considering for instance the normal determined with MDD, it is: N 1 = [0.925, 0.124, -0.355] for interval 1 and: N 2 = [0.872, 0.473, -0.121] for interval 2. Therefore, during the small incursion into the magnetopause which is observed around t = 28s, the normal is different from the normal observed during the large crossing. The two normals are separated by 25°, and the interval 2 normal is closer to the nominal Shue model (which assumes the magnetopause is a paraboloid) than the interval 1 normal.

Nevertheless, one can also observe that, at some points (see t ≈ 22 or t ≈ 29), the results can be significantly different (with fast variations for LNA), while D 1 is not much smaller than unity. A possible reason for these differences may be the use of different current densities: LNA uses the particle current density, whereas MDD is based on the magnetic field. These departures may also indicate that, sometimes, the layer is 1D in the sense of MDD, but not in the sense of LNA. The physical reasons for these discrepancies will be investigated in the next subsection, where the two analysis methods have been tested in a numerical simulation.

Tests of the MDD and LNA methods on a numerical simulation and generalization of MDD

For testing the MDD and LNA methods, we use a 2D numerical PIC simulation published in [START_REF] Dargent | Kinetic simulation of asymmetric magnetic reconnection with cold ions[END_REF]. Note that this simulation of reconnection has no relation with the above experimental case. In this simulation, we have mimicked various spacecraft crossings of the magnetopause layer and treated the data by both the MDD and LNA methods. The crossing used in this paper is shown in Fig. 11 where a map of the magnetic field in the simulation is plotted. The only difference with the real spacecraft data is that the spatial derivatives have been estimated directly from the simulation grid instead of being estimated from the 4-point measurements of the MMS irregular tetrahedron. Fig. 12 shows the results for the crossing shown in Fig. 11, in the same format as Fig. 10, with the same criterion on |∂ t (B)| 2 . It can be seen that MDD determines a normal which is, as expected, close to the y direction, with a clear regular variation which finely fits the shape of the exhaust region in the simulation. It is worth noticing that the B variations are shown to be almost 1D everywhere in the layer, even in the region relatively close to the X point where the field lines are clearly not straight lines.Our LNA result is quite consistent, in general, with this one. Nevertheless, one can once again observe that the two results are not perfectly identical: at some points (see t = 41 -43) where D 1 is very close to unity, the difference between the two results is significant. The LNA result can even include a non negligible z component (not shown), which is inconsistent with the 2D simulation. Although the discrepancies remain generally small, they are to be understood because, for a fully 1D variation, it is clear that j and ∂ t (B) should be strictly tangential and the LNA method should work perfectly. The MDD local normals are plotted also in Fig. 11, where it is clear that the local normal varies along the crossing. These discrepancies point out a weak point in the basic MDD method, which is based on the magnetic field only: when D 1 ≈ 1, it indeed guarantees that the B variations are 1D, so that j is tangential, but it does not guarantee that the other plasma variations are also 1D. In particular, if E variations are not 1D, there is no reason why ∂ t B should be strictly tangential, which is necessary for LNA to work. In low beta regions, one can guess that the magnetic field controls all the other plasma parameters, so that everything is likely to be 1D when the magnetic field is 1D. It is probably the reason why the discrepancies remain quite limited. But in the regions where pressure effects are important (in the central part of the exhaust for instance in reconnection geometries), it is not certain that the 1D variations of B actually ensure the planarity for all the plasma parameters. The fluid equations of momentum, for ions and electrons, clearly show in particular that the variations of the parallel components of the fluid velocities u i and u e are determined by the pressure forces. When these pressure effects are not negligible, the parallel velocities are therefore not constrained by the geometry of the magnetic field variations.

Fortunately, the MDD can easily be generalized. Instead of considering the 3*3 matrix G = ∇B, one can introduce variations of all the needed parameters G = ∇S, where S is a vector of dimension N, including not only the 3 components of B, but also any of the other available parameters: the components of the electric field, those of the ion and electron velocities, those of the pressure tensors, as well as the scalars as the density, etc. In these conditions, G is a 3 * N tensor, but L remains 3*3 and the rest of the method can remain unchanged. A normalization has to be introduced in the computation so that the weight of the different physical quantities is equivalent: the Frobenius norm of ∇B is computed as a function of time, and the magnetic field is normalized by the maximum of the norm over all the interval. And the same is done for the electric field. In the simulation data, such a generalization has been done by just introducing the electric field vector in addition to the magnetic one. The result, which can be compared with the result of Fig. 12 is presented in Fig. 13. One can see that the generalized MDD method allows evidencing a 2D character of the plasma in a small region in the current layer, close to the X point, that was not evidenced by the only B variations. D 1 has more contrasted variations with the non-generalized method, so that the same threshold is now more demanding. This leads to reject some normal determinations in the regions where the discrepancy between the LNA and MDD normals was the most important (with a noticeable z component for the LNA normal in particular) and where D 1 has now smaller values.

Concerning the magnetopause crossing presented in this paper, preliminary tests have been done of the generalization of MDD. They are not presented here because they have not proved yet to be efficient. When applying the same generalization as in the simulation (addition of the E data), the result is not conclusive. The reason seems to be purely experimental: as the calibration of electric antennas is a difficult issue, the precision on the different components of E [START_REF] Ergun | The axial double probe and fields signal processing for the mms mission[END_REF] is not sufficient to calculate safely the tensor ∇E from the four spacecraft measurements: even the basic Maxwell-Faraday law cannot be verified from the data because the differences between spacecraft are dominated by the differences between offsets rather than by the physical differences. The problem is still complicated by the presence, on the magnetospheric side, of very strong electrostatic bursts of short period, which can hardly be eliminated by the smoothing process and which make difficult obtaining the small transverse field induced by the current layers we are interested in. The attempts to use the MDD method modified by [START_REF] Denton | Test of methods to infer the magnetic reconnection geometry from spacecraft data[END_REF] have not allowed hitherto to overcome this difficulty.

Generalizing with the ion velocity V i does not pose similar problems. This has been done, but this test did not lead to conclusive results either: introducing the V i variations does not change significantly the result obtained with B alone. Improving the generalized MDD method to make it efficient with the experimental observations is still a work in progress.

Conclusion and perspectives

For investigating the magnetopause internal structure, one cannot be satisfied with the simplest hypothesis of a perfectly stationary and mono-dimensional layer. We give here evidence of departures from these two simple hypotheses on a magnetopause crossing by MMS. The departure from planarity is particularly investigated, introducing a new single spacecraft method, called LNA, used together with an existing multi-spacecraft method called MDD [START_REF] Shi | Motion of observed structures calculated from multi-point magnetic field measurements: Application to cluster[END_REF]. As LNA can give a reliable result only when the variations are locally 1D, it can indeed be usefully combined with MDD, which allows selecting the intervals where this local 1D hypothesis is verified. We have shown that the basic MDD method, which is based on the B variations only, is not always sufficient for that: even when it indicates variations close to perfectly 1D, the normal provided by LNA can show small but significant differences with the corresponding normal coming from MDD itself. We therefore propose a generalization of MDD using more data. The idea has been tested by adding the E variations to the B ones, with data coming from a numerical simulation: the test has shown that this addition is sufficient for solving, at least partly, the problem. It remains to be investigated more thoroughly with spacecraft data. It is worth emphasizing once again that this paper presents the different methods accessible by MMS for investigating the internal structure of the magnetopause only from a case study: benchmarking these methods and comparing their performances on a statistical basis remain to be done in future studies.

Pending these studies, Table 1 shows that the case presented here is not exceptional and seems rather typical. We analyze six cases in the same way as above, six of them being in the same day as the example of this paper. And we show that the two determinations, LNA and MDD, when restricted to strong criteria for D 1 and for the amplitude of the B variation, are globally consistent, even though they both vary with respect to the "global" MVAB normal (determined in a short interval including the main magnetic gradient). They both show to be often clearly different from this global MVAB determination. The choice of severe criteria has been done here in order to limit as much as possible the effects of non planarity and the role of the superposed turbulence and therefore make the different cases more comparable. However, the results are not perfect in the sense that the distance between the LNA and MDD determinations, which could be expected to be negligible, are generally not smaller than the local variations of each determination, as estimated by the standard deviation of their direction with respect the global MVAB result. This imperfection is likely to be due to the same reason as explained above: using MDD only on the magnetic field does not guarantee the real mono-dimensionality of the physics. Generalizing the method to the electric field should solve this problem if the electric field measurement was accurate enough to allow such a generalization.

The MDD method, contrary to LNA, does not make use of Maxwell equations. In return, it loses the single-spacecraft character of LNA and so part of its locality. is a priori no method that would be strictly single-spacecraft and which would allow to test the local 1D hypothesis with a comparable reliability. Nevertheless, some simplifying hypotheses could be used, in the future, to discard the non-1D regions with some confidence. If one assumes, for instance, that the observed B variations can be approximated locally as stationary in some frame, we must have, in the observation frame:

∂ t (B) = -V • ∇B (7)
where V is the local propagation velocity of the structure. The same property has already been assumed in [START_REF] Shi | Motion of observed structures calculated from multi-point magnetic field measurements: Application to cluster[END_REF], where the propagation velocity of the structures could so be determined. It can be noticed that the red curve plotted in Fig. 5 is an integration of the velocity obtained by this method. The change of slope in the curve around t = 27 indicates a change of the velocity of the boundary and therefore gives a confirmation of the relative back and forth motion of the boundary that was guessed at the beginning of the paper. It seems to also confirm the hypothesis that the flow across the structure is negligible. If true, this may justify Eq. ( 7), the propagation velocity simply being the normal flow velocity.

As soon as the property of Eq. ( 7) is valid, it can easily be shown that the two vectors ∂ t (B) and j are perpendicular to each other when the local variation is 1D, since j = n × ∂ N (B) and ∂ t (B) = -V N ∂ N (B). Checking where the two vectors are perpendicular may provide a test of planarity. This is left for further studies.

As discussed before, the MDD method gives the normal to a one-dimensional boundary, but it can also give information when the problem is 2D. In this case, the eigenvector associated with the largest eigenvalue α 1 does not give much information, but the eigenvector associated with the smallest eigenvalue, α 3 , indicates the direction in which the problem is quasi-invariant. This direction will have to be compared with the direction obtained by other methods such as [START_REF] De Keyser | Empirical reconstruction and long-duration tracking of the magnetospheric boundary in single-and multi-spacecraft contexts[END_REF]. Knowing experimentally the invariant direction may be important for comparing the data with 2D numerical simulations. Of course, α 3 is approximately in the plane perpendicular to N Shue , since the effective normal, given by α 1 is not much different from N Shue . In this plane, investigating the actual direction of α 3 deserves to be explored further. It may provide information, for instance, on the local fluctuations at different scales, whatever their cause: reconnection [START_REF] Aunai | Orientation of the x-line in asymmetric magnetic reconnection[END_REF], Kelvin-Helmholtz [START_REF] Miura | Nonlocal stability analysis of the mhd kelvinhelmholtz instability in a compressible plasma[END_REF], [START_REF] Belmont | Advances in magnetopause kelvin-helmholtz instability studies[END_REF]) or any other phenomenon.

Finally, we have reported in Fig. 3 the intervals where the B variations are mainly 1D (D 1 > 0.98) or 2D (D 3 < 0.05D 2 ) with a colour code. Of course these criteria leave many intervals where the dimension of the problem is not determined, either because the variations are too weak and the concept of dimension is meaningless, either because the dimension of the problem is not close to 1D or 2D. The 2D intervals are concentrated in the region where the spacecraft go back into the magnetopause layer which is reached only in the very small interval around t = 28s. It seems that this incursion is made in a region which is much more complex than the "clean" magnetopause crossing observed at the beginning of the period. The "oscillations" that are seen in the dimension may correspond to the oscillations that are observed on the density. The reason remains to be investigated.

Figure 1 .

 1 Figure 1. Radial distance from the Earth as a function of time: comparison between MMS orbit (blue line) and Shue magnetopause position computed with ACE data.

Figure 2 .

 2 Figure 2. GSE Magnetic field components observed on MMS1, October 16th, 2015, beginning time at 13:05:30.

Figure 3 .

 3 Figure 3. From top to bottom: magnetic field, electron velocity, density, and spectrograms of ions and electrons for the global period studied in the paper. The blue boxes select the regions where the geometry is 1D and the yellow ones the regions where it is 2D (see the discussion at the end of the paper).

Figure 4 .

 4 Figure 4. Hodogram of the magnetic field in the plane tangential to the magnetopause obtained by BV, and its fit. The tangential directions BT1 and BT2 chosen for the plot are those of intermediate and largest variances, but any rotation would not change the interpretation. The axis scales are in nT.

Figure 5 .

 5 Figure 5. Abscissas x(t) along the magnetopause normal, as determined by two different methods (see text). The origin is arbitrary.

Figure 6 .

 6 Figure 6. Comparison of the normal component of the velocity and the density variation (MMS1). The thick lines correspond to the t = 14 -18s interval. The vertical thin lines indicate the limits of the two periods described in the text.

Figure 7 .

 7 Figure 7. Comparison of the main component (B z ) of the magnetic field (left) and computation of the delay (right) between points having the same B z value. The green vertical line is the average delay.

Figure 8 .

 8 Figure 8. The three components of the vector N L N A as determined by LNA without 1D selection in GSE frame, with no test of the significance of the result.

Figure 9 .

 9 Figure 9. The three dimensions resulting of the MDD Analysis as functions of time for the same interval as Fig 8.

Figure 10 .

 10 Figure 10. On top, the D 1 parameter. Below the angle between the normal determined by MDD (in blue) and the reference normal given by the Shue model. The thin lines correspond to D 1 >0.9. The thick lines correspond to D 1 >0.98. In black, the results of the LNA method have been over-plotted for comparison, with the same convention. The intervals selected by dotted lines refer to the text.

Figure 11 .

 11 Figure 11. B z component in the numerical simulation superimposed to the magnetic field lines in the simulation plane. The (x, y) components are those of the 2D simulation box. The straight line indicates the simulated crossing trajectory, with the period of time which is studied below over-lined in green, beginning at the bottom of the simulation box and going in the direction of the increasing y. The small arrows are the MDD local normals determined along the trajectory.

Figure 12 .

 12 Figure 12. Same as Fig. 10 for the crossing in the simulation box shown on Fig. 11. The time is counted from the entrance of the spacecraft in the simulation box which is crossed at constant velocity. The angle is measured with respect to the reference direction, which is here the y direction of the simulation box. The thin lines correspond to D 1 >0.9. The thick lines correspond to D 1 >0.98.

Figure 13 .

 13 Figure 13. Same as Fig. 12 for the crossing in the simulation box shown on Fig. 11 when MDD is replaced by MDD generalized to E field. The three components of the electric field are plotted in the lowest panel for reference.

Table 1 .

 1 There Date θ L N A/MV AB θ M DD/MV AB θ L N A/M DD Comparison of the normals obtained by MDD and LNA on the periods given on the left (the duration is indicated in brackets). The table provides the angles (in degrees) of the two types of normals with respect to MVAB and the angle between them. The statistics are done over all the local normals that satisfy D 1 > 0.99 and ∂ t (B) 2 > 0.5 of its maximum value. The first number corresponds to the mean value and the second one (after ±) corresponds to the standard deviation.

	2015 10 16
	10:20:00 (+120)	20 ± 3	17 ± 8	9 ± 6
	2015 10 16
	10:29:30 (+120)	56 ± 0.5	44 ± 3	12 ± 4
	2015 10 16
	10:36:30 (+120)	33 ± 0.8	21 ± 0.9	12 ± 0.4
	2015 10 16
	10:55:00 (+60)	12 ± 1	11 ± 4	3 ± 1
	2015 10 16
	13:05:30 (+60)	24 ± 2	20 ± 3	7 ± 3
	2017 01 27
	12:05:23 (+70)	35 ± 19	39 ± 14	9 ± 6

Acknowledgments

The authors thank Olivier Le Contel and Laurent Mirioni for their help in dealing with the MMS data and for fruitful discussions. The French involvement on MMS is supported by CNES and CNRS. All the data used are available on the MMS data server: https://lasp.colorado.edu/mms/sdc/public/about/browse-wrapper/


Elena Marchiori

Michèle Sebag
email: sebag@lri.fr

Bayesian learning with local support vector machines for cancer classification with gene expression data

This paper describes a novel method for improving classification of support vector machines (SVM) with recursive feature selection (SVM-RFE) when applied to cancer classification with gene expression data. The method employs pairs of support vectors of a linear SVM-RFE classifier for generating a sequence of new SVM classifiers, called local support classifiers. This sequence is used in two Bayesian learning techniques: as ensemble of classifiers in Optimal Bayes, and as attributes in Naive Bayes. The resulting classifiers are applied to four publically available gene expression datasets from leukemia, ovarian, lymphoma, and colon cancer data, respectively. The results indicate that the proposed approach improves significantly the predictive performance of the baseline SVM classifier, its stability and robustness, with satisfactory results on all datasets. In particular, perfect classification is achieved on the leukemia and ovarian cancer datasets.

Introduction

This paper deals with tumor classification with gene expression data. Microarray technology provides a tool for estimating expression of thousands of genes simultaneously. To this end, DNA arrays are used, consisting of a large number of DNA molecules spotted in a systematic order on a solid substrate. Depending on the size of each DNA spot on the array, DNA arrays are called microarrays when the diameter of DNA spot is less than 250 microns, and macroarrays when the diameter is bigger than 300 microns. DNA microarrays contain thousands of individual DNA sequences printed in a high density array on a glass microscope slide using a robotic instrument. The relative abundance of these spotted DNA sequences in the two DNA and RNA samples may be assessed by monitoring the differential hybridization of the two samples to the sequences on the array. For mRNA samples, the two samples are reverse-transcribed into cDNA, labeled using different fluorescent dyes mixed (red-fluorescent dye Cy5 and green-fluorescent dye Cy3). After these samples are hybridized with the arrayed DNA probes, the slides are imaged using scanner that makes fluorescence measurements for each dye. The log ratio between the two intensities of each dye is used as the gene expression data (cf. [START_REF] Eisen | Dna arrays for analysis of gene expression[END_REF]) expression(gene) = log 2 (int(Cy5)/int(Cy3)), were int(Cy5) and int(Cy3) are the intensities of the two fluorescent dyes.

Four main machine learning tasks are used to analyze DNA microarray data: clustering, e.g. for identifying tumor subtypes, classification, e.g. for tumor diagnostic, feature selection for potential tumor biomarker identification, and gene regulatory network modeling. This paper deals with classification.

Many machine learning techniques have been applied to classify gene expression data, including Fisher linear discriminat analysis [START_REF] Dudoit | Comparison of discrimination methods for the classification of tumors using gene expression data[END_REF], k-nearest neighbour [START_REF] Li | Gene assessment and sample classification for gene expression data using a genetic algorithm/k-nearest neighbor method[END_REF], decision tree, multi-layer perceptron [START_REF] Khan | Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks[END_REF][START_REF] Xu | Artificial neural networks and gene filtering distinguish between global gene expression profiles of barrett's esophagus and esophageal cancer[END_REF], support vector machine (SVM) [START_REF] Brown | Knowledge-based analysis of microarray gene expression data by using support vector machines[END_REF][START_REF] Furey | Support vector machine classification and validation of cancer tissue samples using microarray expression data[END_REF][START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF][START_REF] Noble | Support vector machine applications in computational biology[END_REF], boosting and ensemble methods [START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF][START_REF] Cho | Machine learning in DNA microarray analysis for cancer classification[END_REF][START_REF] Tan | Ensemble machine learning on gene expression data for cancer classification[END_REF][START_REF] Liu | A combinational feature selection and ensemble neural network method for classification of gene expression data[END_REF][START_REF] Dettling | Boosting for tumor classification with gene expression data[END_REF]. A recent comparison of classification and feature selection algorithms applied to tumor classification can be found in [START_REF] Cho | Machine learning in DNA microarray analysis for cancer classification[END_REF][START_REF] Tan | Ensemble machine learning on gene expression data for cancer classification[END_REF].

This paper introduces a method that improves the predictive performance of a linear SVM with Recursive Feature Elimination (SVM-RFE) [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] on four gene expression datasets. The method is motivated by previous work on aggregration of classifiers [START_REF] Breiman | Bagging predictors[END_REF][START_REF] Breiman | Arcing classifiers[END_REF], where it is shown that gains in accuracy can be obtained by aggregrating classifiers built from perturbed versions of the train set, for instance using bootstrapping. Application of aggregration of classifiers to microarray data is described e.g. in [START_REF] Dudoit | Comparison of discrimination methods for the classification of tumors using gene expression data[END_REF][START_REF] Cho | Machine learning in DNA microarray analysis for cancer classification[END_REF][START_REF] Liu | A combinational feature selection and ensemble neural network method for classification of gene expression data[END_REF][START_REF] Dettling | Boosting for tumor classification with gene expression data[END_REF].

In this paper a novel approach is proposed, for generating a sequence of classifiers from the support vectors of a baseline linear SVM-RFE classifier. Each pair of support vectors of the same class are used to generate an element of the sequence, called local support classifier (lsc). Such classifier is obtained by training SVM-RFE on data consisting of the two selected support vectors and all the support vectors of the other class.

The sequence of lsc's provides an approximate description of the data distribution by means of a set of linear decision functions, one for each region of the input space in a small neighbourhoods of two support vectors having equal class label.

We propose to use this sequence of classifiers in Bayesian learning (cf. [START_REF] Mitchell | Machine Learning[END_REF]). The first technique applies Naive Bayes to the transformed data, where an example is mapped into the binary vector of its classification values. The resulting classifier is called Naive Bayes Local Support Classifier (NB-LSC). The second technique applies Optimal Bayes to the sequence of lsc's classifiers. The resulting classifier is called Optimal Bayes Local Support Classifier (OB-LSC).

The two classifiers are applied to four publically available datasets for cancer classification with gene expression. The results show a significant improvement in predictive performance of OB-LSC over the baseline linear SVM-RFE classifier, and a gain in stability. In particular, on the leukemia and ovarian cancer datasets perfect classification is obtained, and on the other datasets performance comparable to the best published results we are aware of.

The rest of the paper is organized as follows. The next two sections describe the baseline and new methods. Sections 4 contains a short description of the data. Section 5 reports results of experiments and discuss them. Finally, the paper ends with conclusive considerations on research issues to be tackled in future work.

Support Vector Machines

This section describes in brief SVM-RFE, the local support classifier construction procedure, and the integration of the resulting classifier sequence in Naive Bayes and Optimal Bayes classification.

SVM

In linear SVM binary classification [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Cristianini | Support Vector machines[END_REF] patterns of two classes are linearly separated by means of a maximum margin hyperplane, that is, the hyperplane that maximizes the sum of the distances between the hyperplane and its closest points of each of the two classes (the margin). When the classes are not linearly separable, a variant of SVM, called soft-margin SVM, is used. This SVM variant penalizes misclassification errors and employs a parameter (the soft-margin constant C) to control the cost of misclassification.

Training a linear SVM classifier amounts to solving the following constrained optimization problem:

min w,b,ξ k 1 2 ||w|| 2 + C m i=1 ξ i s.t. w • x i + b ≥ 1 -ξ i
with one constraint for each training example x i . Usually the dual form of the optimization problem is solved: Maximizing the margin allows one to minimize bounds on generalization error. Because the size of the margin does not depend on the data dimension, SVM are robust with respect to data with high input dimension. However, SVM are sensitive to the presence of (potential) outliers, (cf. [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] for an illustrative example) due to the regularization term for penalizing misclassification (which depends on the choice of C).

min αi 1 2 m i=1 m j=1 α i α j y i y j x i • x j - m i=1 α i such that 0 ≤ α i ≤ C, m i=1 α i y i = 0. SVM requires O(m 2

SVM-RFE

The weights w i provide information about feature relevance, where bigger weight size implies higher feature relevance. In this paper feature x i is scored by means of the absolute value of w i . Other scoring functions based on weight features are possible, like, e.g., w 2 i , which is used in the original SVM-RFE algorithm [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF]. SVM-RFE is an iterative algorithm. Each iteration consists of the following two steps. First feature weights, obtained by training a linear SVM on the training set, are used in a scoring function for ranking features as described above. Next, the feature with minimum rank is removed from the data. In this way, a chain of feature subsets of decreasing size is obtained.

In the original SVM-RFE algorithm one feature is discarded at each iteration. Other choices are suggested in [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], where at each iteration features with rank lower than a user-given theshold are removed. In general, the threshold influences the results of SVM-RFE [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF]. In this paper we use a simple instance of SVM-RFE where the user specifies the number of features to be selected, 70% of the actual number of features are initially removed, and then 50% at each further iteration. These values are chosen after cross-validation applied to the training set.

Local Support Classifiers

We propose to describe the distribution of the two classes by means of a sequence of classifiers, generated from pairs of support vectors ofSVM-RFE. Each of these classifiers, called local support classifier (lsc), is obtained using data generated from two support vectors of the same class, and all support vectors of the other class. In this way, each classifier uses only a local region near the two selected support vectors when separating the two classes. Each classifier generated from two (distinct) support vectors of the same class provides an approximate description of the distribution of the other class given the two selected support vectors.

Before describing the procedure for constructing lsc's, some notation used throughout the paper is introduced.

-D denotes the training set, -c denotes the classifier obtained by training a linear SVM on D, -S p and S n denote the set of positive and negative support vectors of c, respectively, -P air p and P air n denote the set of pairs of distinct elements of S p and S n , respectively.

The following procedure, called LSC, takes as input one (s, s ) in P air p and outputs a linear SVM classifier C s,s by means of the following two steps. An analogous procedure is applied to generate C s,s from pairs (s, s ) in P air n .

When applied to all pairs of support vectors in P air p , P air n , LSC produces a sequence of lsc's. Such sequence of classifiers induces a data transformation, called seq D , which maps example x in the sequence seq D (x) of class values C s,s (x), with (s, s ) in P air p ∪ P air n .

The construction of the sequence oflsc's requires computation that grows quadratically with the number of support vectors. However, this is not a severe problem, since the number of examples, hence of support vectors, is small for this type of data. Furthermore, LSC is applied to each pair of support vectors independently, hence can be executed in parallel.

Naive Bayes and Optimal Bayes Classification

Naive Bayes (NB) is based on the principle of assigning to a new example the most probable target value, given the attribute values of the example. In order to apply directly NB to the original gene expression data, gene values need to be discretized, since NB assumes discrete-valued attributes. Examples transformed using seq D contain binary attributes, hence discretization is not necessary.

Let x be a new example. Suppose seq D (x) = (x 1 , . . . , x N ). First, the prior probabilities p y of the two target values are estimated by means of the frequency of positive and negative examples occurring in the train set D, respectively. Next, for each attribute value x i , the probability P (x i | y) of x i given target value y is estimated as the frequency with which x i occurs as value of i-th attribute among the examples of D with class value y. Finally, the classification of x is computed as the y that maximizes the product

p y N i=1 P (x i | y).
The resulting classifier is denoted by NB-LSC.

Optimal Bayes (OB) classifier is based on the principle of maximizing the probability that a new example is classified correctly, given the available data, classifiers, and prior probabilities over the classifiers.

OB maps example x to the class that maximizes the weighted sum

C s,s w s,s I(C s,s (x) = y),
where w s,s is the accuracy of C s,s over D, and I is the indicator function, which returns 1 if the test contained in its argument is satisfied and 0 otherwise. The resulting classifier is denoted by OB-LSC.

Datasets

There are several microarray datasets from published cancer gene expression studies, including leukemia cancer dataset, colon cancer dataset, lymphoma dataset, breast cancer dataset, NCI60 dataset, ovarian cancer, and prostate dataset. Among them four datasets are used in this paper, available e.g. at http://sdmc.lit.org.sg/GEDatasets/Datasets.html. The first and third dataset contain samples from two variants of the same disease, the second and last dataset consist of tumor and normal samples of the same tissue. Table 1 shows input dimension and class sizes of the datasets. The following short description of the datasets is partly based on [START_REF] Cho | Machine learning in DNA microarray analysis for cancer classification[END_REF].

Leukemia

The Leukemia dataset consists of 72 samples: 25 samples of acute myeloid leukemia (AML) and 47 samples of acute lymphoblastic leukemia (ALL). The source of the gene expression measurements is taken from 63 bone marrow samples and 9 peripheral blood samples. Gene expression levels in these 72 samples are measured using high density oligonucleotide microarrays [START_REF] Ben-Dor | Tissue classification with gene expression profiles[END_REF]. Each sample contains 7129 gene expression levels.

Colon

The Colon dataset consists of 62 samples of colon epithelial cells taken from colon-cancer patients. Each sample contains 2000 gene expression levels. Although the original data consists of 6000 gene expression levels, 4000 out of 6000 were removed based on the confidence in the measured expression levels. 40 of 62 samples are colon cancer samples and the remaining are normal samples. Each sample is taken from tumors and normal healthy parts of the colons of the same patients and measured using high density oligonucleotide arrays [START_REF] Alon | Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays[END_REF].

Lymphoma

B cell diffuse large cell lymphoma (B-DLCL) is a heterogeneous group of tumors, based on significant variations in morphology, clinical presentation, and response to treatment. Gene expression profiling has revealed two distinct tumor subtypes of B-DLCL: germinal center B cell-like DLCL and activated B cell-like DLCL [START_REF] Lossos | Ongoing immunoglobulin somatic mutation in germinal center b cell-like but not in activated b cell-like diffuse large cell lymphomas[END_REF]. Lymphoma dataset consists of 24 samples of germinal center B-like and 23 samples of activated B-like.

Ovarian

Ovarian tissue from 30 patients with cancer and 23 without cancer were analyzed for mRNA expression using glass arrays spotted for 1536 gene clones. Attribute i of patient j is the measure of the mRNA expression of the i-th gene in that tissue sample, relative to control tissue, with a common control employed for all experiments [START_REF] Schummer | Comparative hybridization of an array of 21,500 ovarian cdnas for the discovery of genes overexpressed in ovarian carcinomas[END_REF].

Numerical Experiments

The two classifiers NB-LSC and OB-LSC, described in Section 3.1, are applied to the four gene expression datasets the baseline SVM-RFE algorithm. In all experiments the same value of the SVM parameter C = 10 is used, while the number of selected genes was set to 30 for the lymphoma dataset and 50 for all other datasets. These values are chosen by means of cross-validation applied to the training set.

Because of the small size of the datasets, Leave One Out Cross Validation (LOOCV) is used to estimate the predictive performance of the algorithms [START_REF] Evgeniou | Leave one out error, stability, and generalization of voting combinations of classifiers[END_REF]. Moreover, while the performance of SVM on the Lymphoma dataset is rather scare (possibly due to the fact that we did not scale the data), OB-LSC obtains results competitive to the best results known (see Table 3). Table 3 reports results of OB-LSC and the best result among those contained nine papers on tumor classification and feature selection using different machine learning methods [START_REF] Cho | Machine learning in DNA microarray analysis for cancer classification[END_REF]. Note that results reported in this table have been obtained using different cross-validation methods, mainly by repeated random partitioning the data into train and test set using 70 and 30 % of the data, respectively. Because the resulting estimate of predictive performance may be more biased than the one of LOOCV [START_REF] Evgeniou | Leave one out error, stability, and generalization of voting combinations of classifiers[END_REF], those results give only an indication for comparing the methods. Only the results on the colon dataset from Liu et al 04 and Dettling et al 03 [START_REF] Dettling | Boosting for tumor classification with gene expression data[END_REF][START_REF] Liu | A combinational feature selection and ensemble neural network method for classification of gene expression data[END_REF] are obtained using LOOCV. The methods proposed in these latter papers use boosting and bagging, respectively. The results they obtain seem comparable to OB-LSC.

The results indicate that OB-LSC is competitive with most recent classification techniques for this task, including non-linear methods.

Conclusion

This paper introduced an approach that improves predictive performance and stability of linear SVM for tumor classification with gene expression data on four gene expression datasets.

We conclude with two considerations on research issues still to be addressed. Our approach is at this stage still an heuristic, and needs further experimental and theoretical analysis. In particular, we intend to analyze how performance is related to the number of support vectors chosen to generate lsc's. Moreover, we intend to investigate the use of this approach for feature selection, for instance whether the generated lsc's can be used for ensemble feature ranking [START_REF] Jong | Ensemble feature ranking[END_REF].

  ) storage and O(m 3 ) to solve. The resulting decision function f (x) = w • x + b has weight vector w = m k=1 α k y k x k . Examples x i for which α i > 0 are called support vectors, since they define uniquely the maximum margin hyperplane.

1 .

 1 Let Xp = {s, s }. Assign positive class label to these examples. 2. Let C s,s be the classifier obtained by training a linear SVM on data Xp∪Sn.

Table 1 .

 1 Datasets description

	Name	Tot Positive Negative Genes
	Colon	62	22	40	2000
	Leukemia 72	25	47	7129
	Lymphoma 58	26	32	7129
	Ovarian 54	30	24	1536

Table 2 .

 2 Results of LOOCV: average sensitivity, specificity and accuracy (with standard deviation between brackets).

	Method	Dataset	Sensitivity	Specificity	Accuracy
	SVM-RFE Colon	0.90 (0.3038) 1.00 (0.00)	0.9355 (0.2477)
	NB-LSC		0.75 (0.4385) 1.00 (0.00)	0.8387 (0.3708)
	OB-LSC		0.90 (0.3038) 1.00 (0.00)	0.9355 (0.2477)
	SVM-RFE Leukemia 0.96 (0.20)	1.00 (0.00)	0.9861 (0.1179)
	NB-LSC		1.00 (0.00)	1.00 (0.00)	1.00 (0.00)
	OB-LSC		1.00 (0.00)	1.00 (0.00)	1.00 (0.00)
	SVM-RFE Ovarian	0.7000 (0.4661) 0.9583 (0.2041) 0.8148 (0.3921)
	NB-LSC		1.00 (0.00)	1.00 (0.00)	1.00 (0.00)
	OB-LSC		1.00 (0.00)	1.00 (0.00)	1.00 (0.00)
	SVM-RFE Lymphoma 0.6923 (0.4707) 0.6562 (0.4826) 0.6724 ( 0.4734)
	NB-LSC		1.00 (0.00)	0.6562 (0.4826) 0.8103 (0.3955)
	OB-LSC		1.00 (0.00)	0.8750 (0.3360) 0.9310 (0.2556)

Table 2

 2 reports results of LOOCV. They indicate a statistically significant improvement of OB-LSC over the baseline SVM-RFE classifier, and a gain in stability, indicated by lower standard deviation values. In particular, on the ovarian and leukemia datasets both NB-LSC and OB-LSC achieve perfect classification.

Table 3 .

 3 Comparison of results with best average accuracy reported in previous papers on tumor classification. The type of classifiers considered in the paper are given between brackets. An entry '-' means that the corresponding dataset has not been considered.

	Colon Leukemia Lymphoma Ovarian