Djilali Ameur
email: d.ameur@yahoo.fr

Joanna Dib
email: joannadib2022@yahoo.fr

Soumia Belalia
email: belaliasomia@gmail.com

Isabelle Raspo
email: isabelle.raspo@univ-amu.fr

STATISTICAL STUDY FOR COLLOIDAL SUSPENSIONS COAGULATION

Keywords: Pollution, Polymers, Colloids, Simulation, Smoluchowski equations, Monte Carlo method

Over the years, because of fast growth of population and industrial activities, we have increasingly polluted ou waters. Conventional water treatment processes become indequate with the identification of new contaminants in the water. To overcome these challenges, a particular emphasis is accorded to the implementation of environmentally, ecofriendly and less-energy technologies for water treatment processes. In this context, we investigate wether new chemical treatments employing organic polymers are efficient enough to provide green solutions for waste water recycling. We focus our attention on stabilized colloidal suspensions that show no sedimentation rate. In order to destabilize these colloids, two mechanisms are mainly involved: coagulation and floculation. These mechanisms depend on the aggregation rate between particles modelised by the Smoluchowski's coagulation equation. In this paper, we study the number concentration of coagulated particles and we show more the time step decreases, more the error on moments of order 0 decreases for a certain number of particles.

INTRODUCTION

In order to provide a sustainable development and to avoid the negative impacts on our environment, green chesmistry science and technology offers economically alternatives for chemical water treatments by taking into account the protection of public health and the environment. These techniques are carried out by adding organic polymers which may replace the inorganic coagulants in order to meet the green treatment goals.

The conventional methods of solid-liquid separation such as filtration, sedimentation, centrifugation and flotation, for urban waste process water treatment cannot be used directly on stabilized colloids suspensions that show no sedimentation rate despite the existence of density variation between the two phases. For waste water treatment, difficulties occur when the particle size allows it to remain in suspension in the liquid medium. Classical decanting mechanism is possible when the pretreated water stays in calm waters in the primary settling basin where suspended materials are deposited in the bottom to be scraped and evacuated. But it is unfortunately not enough because the time of decantation of a colloid takes between 2 and 200 years. Colloids are particles of sizes varying between 10 -2 and 1 μm and are generally negatively charged. Actually, colloids have a surface activity which has an appreciable influence on the properties of the aggregate. This phenomena governs the stability of colloidal suspensions which results from the balance between the attractive energies of Van der Waals, which favor aggregation, and the electrostatic repulsion energies. The decantation is carried out by adding polymers whose function is to agglomerate the solid particles in the form of flocs and thus to increase the rate of sedimentation. There are two ways of producing changes in the state of agglomeration of suspensional colloids defined by the terms coagulation and floculation.

In fact, the colloid moves in the the solution with part of its electric double-layer repulsion. The first layer is fixed to the colloid surface and the second is more diffuse. The potential difference between the first layer and the solution is called the zeta potential which governs the displacement of the colloids and their mutual interaction [START_REF] O'brien | A mathematical model for colloidal aggregation[END_REF][START_REF] Henry | Towards a description of particulate fouling: from single-particle deposition to clogging[END_REF]. To destabilize the suspension and promote the agglomeration of the colloids, it is necessary to cross the energy barrier and thus to reduce the zeta potential so that the Van der Waals forces as well as the surface adsorption phenomena become dominant again. This kind of aggregation is known as coagulation. This destabilization is achieved by coagulation polymers of Low Molecular Weight with a high cationic charge in order to neutralize the negative charges of the colloids. Having become true dipoles, these particles are attracted to each other and collide with each other due to the momentum of Brownian motion.

After colloidal suspensions have been destabilized, flocculants polymers are frequently used because of their very High Molecular Weight, i.e. high degree of monomer polymerization. These flocculants are extremely effective in bonding the micro-flocs formed during coagulation by using long polymer chains. The polymer chain in solution is partially adsorbed on a particle beyond its double layer and when other particles are close enough, the elongate chain is adsorbed on its surface and creates bridging between the two particles. In general, the efficiency of the polymer depends on the molecular chain lengh, because more the chain is longer more the probability of creating bridging is greater.

MATHEMATICAL MODEL

For a given polymer, bridging mainly depends on two parameters: the number of sites available for adsorption on the colloids surface and the rate of polymers agglomeration. One can consider that the probability of two polymer chains fusion is proportional to their lengths sum.

The beginning of the classical understanding of colloidal aggregation stem from the work of Smoluchowski on coagulation processes in colloids [START_REF] Smoluchowski | Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lӧsungen[END_REF], a population equilibrium equation which describes in statistical chemical physics the agglomeration which follows from the assumption that the collisions are binary and that fluctuations in density are sufficiently small so that the collisions occur at random. The Smoluchowski equation is of universal use in the fields of colloid chemistry, aerosol dynamics and Nanotechnology Science. It is to say that the collision rate per unit length is due to the Brownian motion and the particle agglomeration rate depends on the particles contact probability and the effectiveness of the latter. Numerical resolution of the Smoluchowski equation is obtained employing the Monte Carlo technique [START_REF] Ameur | Modélisation analytique et simulation numérique par la méthode de Monte Carlo d'un écoulement de gaz dans des micro-canaux[END_REF] and its alternatives methods.

If the polymers are assimilated to their chain length and if we consider that two polymers of respective lengths i and j have a probability K(i,j) to coagulate to give a particle of length i+j, a balance of continuous equation makes it possible to obtain the integro-differential equation [START_REF] Babovsky | On a Monte Carlo scheme for Smoluchowski's coagulation equation[END_REF]: where c(i,t) denotes the number of polymers of length i at time t and K(i,j) the coagulation kernel. Physically, we assume that the number and total length of the polymers are finite and that K(i,j) is nonnegative and symmetric.

The first term on the right-hand side represents the rate of agglomeration of polymers of length k by agglomeration of any polymers pair of length i and j, such that i+j =k. The second term accounts for the loss of polymers of length k by agglomeration with any other polymers chains.

Rather than approximating the density, we shall approximate the mass density which satisfises the following equation: An analytical solution of the problem is given by: 𝑐(𝑖, 𝑡) = 𝑒 -𝑡 𝑖 𝑖-1 𝑖! (1 -𝑒 -𝑡 ) 𝑖-1 𝑒 -𝑖(1-𝑒 -𝑡 ) , 𝑖 ∈ ℵ * , 𝑡 > 0.

𝜕𝑐∂ 𝑡 (𝑖, 𝑡) = 1 2 ∑ 𝐾 ̌(𝑖 -
We compute the solution up to time T=2 for different time steps varying from ∆t = T 100 to ∆t = T 1000 and for several particle numbers varying from N=10 4 to N=10 5 and we compare the approximate moment of order 0 which corresponds to the number concentration of coagulated particles which decreases with the growth of the particle size:

𝑀 0,𝑁 (𝑡) = ∑ 1 𝑖 𝑁,𝑛 (𝑘) 𝑁 𝑘=1
to the exact moment of order 0 over the interval [0, 2]:

𝑀 𝑝 (𝑡) = ∑ 𝑖 𝑝 ∞ 𝑖=1 𝑐(𝑖, 𝑡) = ∑ 𝑖 𝑝-1 ∞ 𝑖=1 𝑐(𝑖, 𝑡) 𝑤ℎ𝑒𝑟𝑒 𝑀 0 (𝑡) = 𝑒 -𝑡 .

SIMULATION RESULTS

There has recently been a great deal of interest in the computer simulation of aggregation phenomena in order to confirm the theoratical studies. In [START_REF] Dib | Aggregation Rate Simulation of Water Soluble Polymers for Urban Waste Water Process Treatment[END_REF], we have studied the moments of order 0 and 2 and we have shown that the method sems to be more accurate approximating numerically the first moment than the second one usually used as an index to characterize the total light scattered, which increases with the growth of particle size and polydisperity. In this paper, we focus our study on the moment of order 0 and we show that the time step ∆t influence the error on moment in such a manner that the error varies between 10 -2 and 10 -4 for ∆t varying between 1 100 and 1 300 . Most of the results reported in this paper were obtained from simulations in which the effective total number of colloidal suspensions was varying between 10 4 and 10 5 .

Simulations were performed for a time step Δt varying between 0,000612 0,000632 0,000356 0.000890 0.000073

Table 1: Variation of the error on moment of order 0 in fonction of the number of particles and the time step.

The following simulations show the validity of the Smoluchowski theory approach with comparing results of computer simulations.

  𝑖 -𝑗, 𝑗) 𝑐(𝑖 -𝑗, 𝑡) 𝑐(𝑗, 𝑡) , 𝑗) = 𝑖 + 𝑗 and 𝑐(𝑖, 0) = 𝑐 𝑖 (0) = { 1 𝑖𝑓 𝑖 = 1, 0 𝑒𝑙𝑠𝑒.

1 𝑗=1-

 1 𝑗, 𝑗) 𝑐(𝑖 -𝑗, 𝑡) 𝑐(𝑗, 𝑡) 𝑖-∑ 𝐾 ̌(𝑖, 𝑗) 𝑐(𝑖, 𝑡) 𝑐(𝑗, 𝑡), ∞ 𝑗=1 ∀𝑖 ∈ ℵ * , 𝑡 > 0 , 𝑐(𝑖, 0) = 𝑐̌𝑖(0), where 𝑐(𝑖, 𝑡) = 𝑖 𝑐(𝑖, 𝑡) and 𝐾 ̌(𝑖, 𝑗) = 𝐾(𝑖, 𝑗) 𝑗 .

Figure 1 :

 1 Figure 1: Comparison between analytical and numerical moments of order 0 (left) and the error on moments of order 0 (right) corresponding to N=7x10 4 for ∆𝑡 = T 300 .

Figure 2 :

 2 Figure 2: Comparison between analytical and numerical moments of order 0 (left) and the error on moments of order 0 (right) corresponding to N=8x10 4 for ∆𝑡 = T 400.

Figure 3 :Table 2 :

 32 Figure 3: Comparison between analytical and numerical moments of order 0 (left) and the error on moments of order 0 (right) corresponding to N=10 5 for ∆𝑡 = T 500 . Other simulations were carried out for various time steps varying between ∆𝑡 = 1 10 and ∆𝑡 = 1 60 corresponding to a number of particles of order 10 4 , in order to show the influence of the time step on the error on moment of order 0. One can see that more the time step decreases, more the error decreases solution and the method is validated. ∆𝑡 Error on 𝑴 𝟎,𝑵 ∆𝑡 Error on 𝑴 𝟎,𝑵

  The following table shows that considering a certain number of particles N, more the time step ∆𝑡 decreases, more the error on moments of order 0 𝑴 𝟎,𝑵 decreases. These results prove the precision of the method employed in this paper.

					1 100	and	1 300	and five effective total number of
	particles. Error on 𝑴 𝟎,𝑵	N=3x10 4	N=5x10 4	N=7x10 4	N=8x10 4	N=10 5
	∆𝑡 =	1 100	0,013736	0,004121	0,002829	0,003407	0.003148
	∆𝑡 =	1 200	0,002576	0,002876	0,000602	0,002224	0.001073
	∆𝑡 =	1 300

CONCLUSIONS

In this paper, we have analysed an approach for solving Smoluchowski's coagulation equation using the Monte Carlo method in order to simulate the aggregation of polymers employed in the waste water process treatment.

Our study provides numerical evidence that the numerical method seems to be accurate approximating the zeroth moment representing the particle number concentration, which decreases with the growth of the particle size. We have shown that the time step influence the error on moment in such a manner that considering a certain number of particles, more the time step decreases, more the error on moments of order 0 decreases.


Smahi Housseyne
email: housseynsmahi@gmail.com

Djilali Ameur
email: d.ameur@yahoo.fr

Joanna Dib
email: joannadib2022@yahoo.fr

Isabelle Raspo
email: isabelle.raspo@univ-amu.fr

SIMULATION OF CONVECTIVE MASS TRANSFER IN A NEAR-CRITICAL BINARY MIXTURE ABOVE AN ADSORBENT BOUNDARY

Keywords: Mass transfer, Adsorption, Supercritical carbon dioxide (scCO2), mixture fluid, Piston effect GLOSSARY

In this work, we present a study of the adsorption for dilute solute (Naphthalene) in supercritical CO2 enclosed in a rectangular cavity. The effect of proximity to the critical point on mass transfer is numerically investigated. The numerical model is based on the resolution of the Navier-Stokes equations coupled with the energy and mass diffusion balances. The velocity field, the mass fraction perturbation as well as the temporal evolution of the mean Sherwood number is presented for various temperatures. The results show that close to the mixture critical point, the homogeneous adsorption of the solute is obtained. Additional key findings include that closer we get to the critical point, the better is the mass transfer on the adsorbent plate.

INTRODUCTION

The supercritical fluid (SCF) is a fluid for which temperature and pressure are larger than those of the critical point. In this zone of the phase diagram, the physical properties of the fluid (density, viscosity, diffusivity) are intermediate between those of liquids and gases [START_REF] Ameur | Numerical simulation of the Poiseuille-Rayleigh-Bénard instability for a supercritical fluid in a mini-channel[END_REF][START_REF] Ameur | Thermoconvective instabilities of 2D Poiseuilles-Rayleigh-Benard for supercritical fluids in micro/macro-channels[END_REF]. Supercritical fluid deposition (SFD) is a promising technique aiming to prepare supported metallic nanoparticles or metallic films and has been the subject of several researchs and publications [START_REF] Wannassi | Numerical study of non-isothermal adsorption of Naphthalenein supercritical CO2: Behavior near critical point[END_REF][START_REF] Erkey | Preparation of metallic supported nanoparticles and films using supercritical fluid deposition[END_REF][START_REF] Bozbag | Supercritical deposition: Current status and perspectivesfor the preparation of supported metal nanostructures[END_REF]. In this technique, supercritical CO2 (scCO2) has been used in most of the depositions because it is non-flammable and non-toxic. Using supercritical fluids (SCFs) to deposit metals and metaloxides on surfaces has been attracting due to the attractive properties of supercritical fluids. Supercritical fluid deposition (SFD) in particular using supercritical carbon dioxide (scCO2) is attractive because no liquid the waste is generated, no solvent residue is left on the substrate and the mass transfer rates are fast compared to the liquids [START_REF] Erkey | Preparation of metallic supported nanoparticles and films using supercritical fluid deposition[END_REF]. In addition, the rate and the velocity of the substrate deposition can be controled by the changing of the pressure and temperature of solvent presented here by the supercritical carbon dioxide (scCO2). In this paper, Naphthalene was chosen as a model solute because its phase equilibria with CO2 has been thoroughly studied [START_REF] Wannassi | Numerical study of non-isothermal adsorption of Naphthalenein supercritical CO2: Behavior near critical point[END_REF][START_REF] Raspo | Heterogeneous reaction induced by the piston effect in supercritical binary mixtures[END_REF][START_REF] Raspo | Fast mass transfer at a solid-supercritical fluid interface by piston effect[END_REF]. There are extensive data available for this system that have been confirmed.

FIGURE 1 Physical configuration

The physical model consists of a dilute solute (Naphthalene-scCO2) mixture enclosed in a rectangular cavity with a heated and adsorbent plate on the bottom boundary so showed in figure 1 of height H=1mm. An adiabatic boundary condition was applied to the non-reactive vertical and upper walls. Its temperatures are maintained fixed to Ti. The fluid is considered in thermodynamic equilibrium at a constant temperature Ti slightly above a particular mixture critical temperature Tcm= 307.65 K and Ti=(1+ε), where ε=( Ti-Tcm)/ Tcm defines the dimensionless proximity to the critical point. The critical propertie Tcm corresponds to the Lower Critical EndPoint of the mixture and is slightly above the critical point of CO2. A weak gradually heating is then applied at the solid bottom plate (y = 0 for H<x<L-H).

MATHEMATICAL MODEL

The mathematical model is based on the 2D time-dependent and compressible Navier-Stokes equations, coupled with energy and mass diffusion equations including the supplemental Peng-Robinson equation of state. In order to reduce computational costs, these equations are solved in the framework of the low Mach number approximation [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF]. Thus, the total pressure is split into two parts: a homogeneous thermodynamic part Pth(t), which appears in the equation of state and in the energy equation and only depends on time t, and a non-homogeneous dynamic part Pdyn(x, y, t), appearing in the momentum equation and which varies with time and space. The dimensionless governing equations are therefore:

𝜕𝜌 𝜕𝑡 + 𝜕 𝜕𝑥 (𝜌𝑉) = 0 (1)
𝜌 𝜕𝑉 𝜕𝑡 + 𝜌𝑉. ∇𝑉 = -∇𝑃 𝑑𝑦𝑛 + 1 𝑅𝑒 ∆𝑉 + 1 3𝑅𝑒 ∇(∇. 𝑉) + 1 𝐹𝑟

𝜌

(2)

𝜌 𝜕𝑇 𝜕𝑡 + 𝜌𝑉. ∆𝑇 = - 𝐶 𝑉0 𝐶 𝑉𝑖 (g 0 -1) [𝑃 𝑡ℎ -𝑇 ( 𝜕𝑃 𝑡ℎ 𝜕𝑇 ) 𝜌,𝜔 ] (∇. 𝑉) + g 𝑅𝑒𝑃𝑟 ∇. (l * ∇𝑇) -[(𝑈 ̅ 2 * -𝑈 ̅ 1 * ) + 𝐶 𝑉0 𝐶 𝑉𝑖 (g 0 -1) × (𝑃 𝑡ℎ -𝑇 ( 𝜕𝑃 𝑡ℎ 𝜕𝑇 )) 𝜌,𝜔 (𝑉 ̅ 2 * -𝑉 ̅ 1 * )] × 1 (g-1) 2 𝐿𝑒𝜃(𝜔) ∇. (𝜌𝐷 21 * ∇𝜔) (3)
𝜌 𝜕𝜔 𝜕𝑡 + 𝜌𝑉. ∇𝜔 = 1 (g-1) 2 𝐿𝑒 ∇( 𝜌𝐷 21 * ∇𝜔) (4)
V is the velocity of components u and v in the x-and y-directions respectively, w is the mass fraction, g is the ratio of the isobaric and isochoric specific heats calculated from the equation of state with g 0 and Cv0

corresponding to the values for a perfect gas (g 0 = 1. ) 𝑇 𝑐𝑚 is the sound speed and VPE=H/ti is the characteristic velocity of the piston effect [START_REF] Raspo | Heterogeneous reaction induced by the piston effect in supercritical binary mixtures[END_REF][START_REF] Raspo | Fast mass transfer at a solid-supercritical fluid interface by piston effect[END_REF][START_REF] Accary | An adaptation of the low Mach number approximation for supercritical fluid buoyant flows[END_REF]. It must be noted that the divergence of the infinite-dilution properties of the solute is not specific of Naphthalene but it is a universal behavior for dilute mixtures [START_REF] Wannassi | Numerical study of non-isothermal adsorption of Naphthalenein supercritical CO2: Behavior near critical point[END_REF].

RESULTS

Simulations were performed for a heating ΔT=10mK and three initial temperatures Ti=308.15K, 309.15K and 318.15K, which correspond to distances to the mixture critical point equal to 0.5K, 1.5K and 10.5K respectively. The mixture is saturated, which means that the initial mass fraction wi corresponds to the solubility of Naphthalene in CO2. The values of the various properties and transport coefficients at the initial states considered are given in Table 1. As already said in section 2, since the binary mixture is very dilute because of the low solubility of the solute, we can assume that the transport coefficients of the mixture are those of the pure solvent. Thus, the thermal expansion coefficient β and the ratio of heat capacities γ were calculated from the Peng-Robinson equation for pure CO2. 4

The Damköhler number was fixed to Da=10 -5 , which is a pertinent value for an adsorption reaction. Figure 2 shows the isolines of temperature and mass fraction perturbation at two times for the three initial temperatures. The first time corresponds to 2.5tPE for all the cases. The solutions obtained for the various values of Ti are similar. The temperature field at this time is governed by the Piston effect, as in pure CO2 [START_REF] Amiroudine | Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below[END_REF], with three distinct zones: a hot boundary layer along the heated plate, cold boundary layers along isothermal walls and the cavity bulk which is homogeneously heated by the Piston effect. In addition, as in pure CO2, the thermal boundary layers are thinner when approaching the critical point because of the strong decrease of thermal diffusion.

Ti=308.15K

Ti=309.15K Ti=318.15K In the same way, the isolines of the mass fraction perturbation reveal the existence, along the adsorbent plate, of a thin boundary layer which becomes thinner when approaching the critical point because of the decrease of the mass diffusion coefficient D21 (Table 1). Figure 2 also presents the steady solution obtained for each initial temperature. Concerning the mass fraction, it must be underlined that, although the isoline map does not evolve anymore, the values of w-wi keeps on increasing while the adsorption reaction still goes on. As expected, the temperature field is very different depending on the proximity to the critical point. Indeed, because of the diverging behavior of the thermal expansion coefficient β (Table 1), the Rayleigh number strongly increases as approaching the critical temperature: Ra=3683 for Ti=318.15K, Ra=21938 for Ti=309.15K and Ra=34139 for Ti=308.15K. As a consequence, the convective flow is much more intense for Ti=308.15K (Fig. 3.(a)), allowing a better mixing inside the cavity. The distribution of the convection cells directly affects the mass fraction field as it is revealed in figure 3 showing the velocity field together with the isolines of w-wi. The jet flow between contrarotative rolls is much weaker for Ti=318.15K (Fig. 3.(b)) allowing the development of a high mass fraction "plume" rising at the middle of the heated plate (Fig. 2). For Ti=308.15K, in contrast, the very intense convection flow restricts the diffusion of the solute in the boundary layer along the plate. The mass fraction field obtained for Ti=309.15K is intermediate between those previously described (Fig. 2). The convection flow has also a strong influence on the adsorbed amount of solute. Figure 4 presents the profile of the relative mass fraction perturbation (w-wi)/wi on the bottom boundary for the three initial temperatures.

Although the maximum amount adsorbed on the plate increases when moving away from the critical point, the profiles reveal that the adsorption of the solute is strongly inhomogeneous for Ti=318.15K, with a peak of adsorption in the middle of the plate. On the other hand, close to the critical point (Ti=308.15K), the peak at x=L/2 is less pronounced, leading to a more homogeneous adsorption of the solute and, therefore, a more homogeneous material in the framework of SFD processes for example.

Finally, the effect of the proximity to the critical point on the mass transfer at the adsorbent plate was considered. For this purpose, we calculated the mean Sherwood number on the plate. This dimensionless number, which characterizes the mass transfer at a boundary, is defined by: Figure 5 shows the temporal evolution of the Sherwood number for the three initial temperatures. A Strong value of Sh is observed at the very beginning of the simulations because of the increase of the plate temperature from Ti to Ti+ΔT. Then the Sherwood number rapidly decreases to reach its steady state value. The drop observed between around 5s and 10s is due to the impact of a large cold thermal plume on the heated plate. Figure 5 reveals that the Sherwood number in the steady solution strongly increases when approaching the critical point, with a value of Sh about 1.5 times larger at Ti=308.15K than the one at Ti=318.15K. In other words, getting closer to the mixture critical point enables a better mass transfer at the adsorbent plate. FIGURE5. Temporal evolution of the mean Sherwood number for the three initial temperatures.

CONCLUSIONS

In this paper, we presented a numerical study of heat and mass transfer in a dilute near-critical binary mixture enclosed in a rectangular cavity with a heated and adsorbent plate on the bottom boundary. The Naphthalene-CO2 mixture was chosen as a model mixture because of the substantial experimental data on solubility of Naphthalene, allowing a correct determination of binary parameters in the mixing rules of the equation of state. However, since the diverging behaviors of transport coefficients near critical points are universal, we believe that the conclusions of our study could be applied to any binary mixture exhibiting the same kind of phase diagram. The results showed that, because of the divergence of the thermal expansion coefficient, a more homogeneous adsorption profile can be obtained for a same value of the heating when moving closer to the mixture critical point. In addition, approaching the critical point also leads to a better mass transfer at the adsorbent plate. These two phenomena can be of great interest for SFD processes.

FIGURE 2 .

 2 FIGURE 2. Isolines of temperature (left column) and mass fraction perturbation w-wi (right column) for the three initial temperatures.

FIGURE 3 .

 3 FIGURE 3. Velocity field together with the isolines of w-wi in the central part of the cavity (above the heated plate) for (a) Ti=308.15K and (b) Ti=318.15K.

  ) 𝑖 (𝑤(𝑥, 0, 𝑡) -𝑤(𝑥,

  FIGURE4. Profiles of the relative mass fraction perturbation on the bottom boundary for the three initial temperatures.

TABLE 1 .

 1 Initial values of temperature, mass fraction and transport coefficients.

	Ti	wi	i	λi	CVi	(D21)i	βi	γ
	308.15 7.6751×10 -3 3.3402×10 -5 9.6172×10 -2 1306.27 2.1969×10 -8 2.3134×10 -1 16.75
	309.15 8.0713×10 -3 3.3448×10 -5 9.1322×10 -2 1269.32 2.2009×10 -8 1.7465×10 -1 13.95
	318.15 1.2557×10 -2 3.3870×10 -5 6.6804×10 -2 1074.57 2.2368×10 -8 5.4955×10 -2 6.51


Catherine Joly

Christophe Bliard

Eliane Espuche

M Escoubes

Acetylated

Acetylated Wheat Starch Membranes: Permeability Measurements

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


est destinée au dépôt

On the instability of radially sheared and axially stratified vortex flows

The stability of rotating stratified vortices has been the subject of a large amount of studies as it concerns fluid mechanics in general but geo-and astophysical fluid dynamics in particular. Simple inviscid stability criteria have been obtained for homogeneous axisymmetric rotating vortices characterized by their velocity profile. These criteria address centrifugal but also shear instability. Adding linear stratification with respect to the rotation axis does not affect these criteria but gives rise to a new type of instability, the so called stratorotational instability. Using a thermally stratified Taylor-Couette experiment and two numerical codes applying either the Boussinesq (BA) or the low Mach number approximation (LMA) we study this instability for different azimuthal flow profiles (ranging from quasigalactic, Keplerian to Taylor-Couette type flows) and axial temperature gradients. The new experiments confirm earlier theoretical findings but also show interesting deviations due to the fact that the small gap approximation does not hold for the used experimental design. We also address the question whether for axial temperature differences of about 10K and larger the Boussinesq approximation still holds. For this purpose we compare solutions performed with the BA and the LMA. In the future we plan to do more detailed numerical simulations using massively parallel computing.

  On the instability of radially sheared and axially stratified vortex flows Gabriel Meletti de Oliveira (1), Isabelle Raspo (2), Torsten Seelig (1), Stephane Viazzo (2), Anthony Randriamampianina (2), Stephane Abide (3), Andreas Krebs (1), and Uwe Harlander (1)(1) Brandenburg University of Technology (BTU) Cottbus, Aerodynamics and Fluid Mechanics, Cottbus, Germany (haruwe@b-tu.de), (2) Laboratoire M2P2, UMR 7340 CNRS, AMU, Centrale Marseille, France, (3) LAMPS Laboratory of Mathematics and Physics, Université de Perpignan, France

Gabriel Meletti de Oliveira, Isabelle Raspo, Torsten Seelig, Stéphane Viazzo, Anthony Randriamampianina, Stéphane Abide, Andreas Krebs, Uwe Harlander To cite this version: Gabriel Meletti de Oliveira, Isabelle Raspo, Torsten Seelig, Stéphane Viazzo, Anthony Randriamampianina, et al.. On the instability of radially sheared and axially stratified vortex flows. European Geosciences Union General Assembly 2018, Apr 2018, Vienne, Austria. hal-02117752


Guillaume Damiand
email: damiand@sic.univ-poitiers.fr

Martine Dexet
email: dexet@sic.univ-poitiers.fr

Pascal Lienhardt
email: lienhardt@sic.univ-poitiers.fr

Eric Andres
email: andres@sic.univ-poitiers.fr

Martine Dexet-Guiard

Removal and contraction operations to define combinatorial pyramids: application to the design of a spatial modeler

Keywords: removal, contraction, irregular pyramids, n-G-map, modeler, discrete reconstruction

Introduction

Hierarchical representations are the bases of several applications in the field of discrete imagery. Our first goal is the study of basic problems related to the definition of hierarchical structures. We study the definition of removal and contraction operations, from which hierarchical structures can be derived (cf. [START_REF] Kropatsch | Abstraction pyramids on discrete representations[END_REF]).

Then we study a first experimentation of these notions. A four level hierarchical structure is proposed in order to handle coherent discrete and continuous representations of geometrical objects. The other levels of the structure correspond to other representations which make a link between the previous ones, mainly a discrete analytical representation.

Many works deal with regular (cf. e.g. [START_REF] Burt | Segmentation and estimation of image region properties through cooperative hierarchical computation[END_REF]) or irregular (cf. e.g. [START_REF] Meer | Stochastic image pyramids[END_REF][START_REF] Montanvert | Hierarchical image analysis using irregular tesselations[END_REF][START_REF] Jolion | The adaptive pyramid : a framework for 2d image analysis[END_REF]) image pyramids for multi-level image analysis and processing (cf. also [START_REF] Kropatsch | Abstraction pyramids on discrete representations[END_REF]). For irregular pyramids, it is necessary to handle a (topological) representation and basic operations, for instance dual graphs and removal and contraction operations for 2D images [START_REF] Kropatsch | Building irregular pyramids by dual graph contraction[END_REF][START_REF] Kropatsch | Building irregular pyramids by dual-graph contraction[END_REF]. Similar problems about multi-level representations arise also in geometric modeling (e.g. for CAD applications, architectural or geological modeling,. . . [START_REF] Bonneau | Hierarchical decomposition of datasets on irregular surface meshes[END_REF][START_REF] De Floriani | A hierarchical boundary model for solid object representation[END_REF][START_REF] Aubry | Building hierarchical solid models from sensor data[END_REF]).

Our goal is to build a theoretical framework for the definition and handling of n-dimensional irregular pyramids, in order to get coherent definitions of data structures and operations for any dimension. So we study the definition of removal and contraction of i -dimensional cells within n-dimensional objects (see also [START_REF] Damiand | Removal and contraction for n-dimensional generalized maps[END_REF]), in order to rigorously define the relations between two consecutive levels of a pyramid.

We choose to study the definitions of removal and contraction operations for n-dimensional generalized maps, since this notion enables us to unambiguously represent the topology of quasi-manifolds, which is a well-defined class of subdivisions [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasimanifolds[END_REF]. Note that several models based on combinatorial maps [START_REF] Edmonds | A combinatorial representation for polyhedral surfaces[END_REF][START_REF] Jacques | Constellations et graphes topologiques[END_REF][START_REF] Cori | Un code pour les graphes planaires et ses applications[END_REF] have been proposed for handling 2-dimensional [START_REF] Brun | Segmentation d'images couleur à base topologique[END_REF][START_REF] Fiorio | Approche interpixel en analyse d'images : une topologie et des algorithmes de segmentation[END_REF] and 3-dimensional segmented or multi-level images [START_REF] Braquelaire | A topological structuring for aggregates of 3d discrete objects[END_REF][START_REF] Bertrand | Topological encoding of 3d segmented images[END_REF][START_REF] Bertrand | Topological map: minimal encoding of 3d segmented images[END_REF][START_REF] Braquelaire | 3d split and merge with 3maps[END_REF][START_REF] Damiand | Définition et étude d'un modèle topologique minimal de représentation d'images 2d et 3d[END_REF]. For instance, Brun and Kropatsch revisit works about graphs [START_REF] Brun | Dual contraction of combinatorial maps[END_REF][START_REF] Brun | Pyramids with combinatorial maps[END_REF][START_REF] Brun | The construction of pyramids with combinatorial maps[END_REF][START_REF] Brun | Contraction kernels and combinatorial maps[END_REF] in order to define 2-dimensional combinatorial map pyramids.

We prefer to use generalized maps instead of combinatorial maps, since their algebraic definition is homogeneous and so we can provide simpler definitions of data structures and operations with generalized maps, and so more efficiency for the conception of softwares (Note that several kernels of geometric modeling softwares are based upon data structures derived from this notion). Last, we know how to deduce combinatorial maps from generalized maps, so the results presented in this paper can be extended for combinatorial maps. Precise relations between generalized and combinatorial maps, and other classical data structures are presented in [START_REF] Lienhardt | Topological models for boundary representation: a comparison with n-dimensional generalized maps[END_REF].

We propose a first experiment of the definition of n-dimensional irregular pyramids to define a four level hierarchical structure. This pyramid is the kernel of a discrete geometrical modeler. The goal of this modeler is to handle both discrete and continuous representations of geometrical objects. Operations in the discrete or in the continuous world are quite different. Many of these oper-ations are best performed in one or the other world (morphological operations, intersections, transformations, . . . ). Our goal is to perform an operation is the world where it is best defined. There are applications where both worlds can be needed. One that comes to mind is medical imaging where acquired data is discrete. However, the popularity of techniques such as the Marching cubes show that there many applications in medical imaging where a continuous representation of discrete data is useful.

The particularity of our modeler is that for each object a discrete and a Euclidean representation coexist in the modeler. These representations must be consistent: the discretization1 of the Euclidean representation must be equal to the discrete object. Thus, we also want to be able to update all representations in case of modification in order to ensure the coherence. For instance, each operation applied on the discrete representation (such as Boolean or morphological operations) must automatically update the Euclidean representation and vice versa.

We add some intermediate levels between the discrete and the Euclidean representations in order to divide the updating operation into several basic steps. So, we define a pyramid in order to represent all different representations in a same data structure. This is a four level pyramid: level 0 is the discrete representation (pixels); level 1 is the representation of the borders of regions; level 2 is the discrete analytical representation and level 3 is the Euclidean representation.

We recall in Sec. 2 the notion of generalized maps. Then we define removal of one i -dimensional cell in Sec. 3, and contraction by duality (in Sec. 4). In Sec. 5, we show that it is possible to simultaneously remove and contract several cells of any dimensions, if and only if cells are disjoined. We present in Sec. 6 the four levels pyramid, basis of the discrete geometrical modeler. Then, we present all conversion operations between consecutive levels of the pyramid in Sec. 7. Last we conclude and give some perspectives in Sec. 8.

Generalized Maps Recall

An n-dimensional generalized map is a set of abstract elements, called darts, and applications defined on these darts: (1) B is a finite set of darts;

Definition 1 (Generalized map) Let n ≥ 0. A n-dimensional generalized map (or n-G-map) is G = (B, α 0 , . . . , α n ) where:
(2) ∀i, 0 ≤ i ≤ n, α i is an involution2 on B;

(3) ∀i, j, 0 ≤ i < i + 2 ≤ j ≤ n, α i α j is an involution.

Let G be an n-G-map, and S be the corresponding subdivision. Intuitively, a dart of G corresponds to an (n + 1)-tuple of cells (c 0 , . . . , c n ), where c i is an i -dimensional cell that belongs to the boundary of c i+1 (cf. [START_REF] Brisson | Representing geometric structures in d dimensions: topology and order[END_REF] and Fig. 1). α i associates darts corresponding with (c 0 , . . . , c n ) and (c ′ 0 , . . . , c ′ n ), where c j = c ′ j for j = i, and c i = c ′ i (α i swaps the two i -cells that are incident to the same (i -1) and (i + 1)-cells). When two darts b 1 and b 2 are such that b

1 α i = b 2 (0 ≤ i ≤ n), b 1 is said i-sewn with b 2 .
G-maps represent cells in an implicit way:

Definition 2 (i-cell) Let G be an n-G-
′ 0 = 4(α 0 α 1 )α 0 = 1 ∈ B S . <> N -{i} (b) =< α 0 , . . . , α i-1 , α i+1 , . . . , α n > (b)
Intuitively, an i -cell is the set of all darts which can be reached starting from b, by using any combination of all involutions except α i . The set of i -cells is a partition of the darts of the G-map, for each i between 0 and n. Two cells are disjoined if their intersection is empty, i.e. when no dart is shared by the cells. More precisions about G-maps are provided in [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasimanifolds[END_REF].

Removal

Intuitively and in a general way for an n-dimensional space, the removal of an i -cell consists in removing this cell and in merging its two incidents (i+1)-cells: so removal can be defined for 0 . . . (n -1)-cells.

Dimension 1: 0-Removal

For dimension 1, only the 0-removal exists, which consists in removing a vertex and in merging its two incident edges. Let C =< α 1 > (b) be a vertex, let Cα 0 be the "neighbor" darts of C for α 0 , i.e. Cα 0 = {b ′′ | ∃b ′ ∈ C such that b ′ α 0 = b ′′ }, and let B S = Cα 0 -C be the "neighbor" darts of C for α 0 that do not belong to C (see Fig. 2). The G-map resulting from the 0-removal of C is obtained by redefining α 0 for the darts of B S as follows:

∀b ′ ∈ B S , b ′ α ′ 0 = b ′ (α 0 α 1 ) k α 0 ,
where k is the smallest integer such that b ′ (α 0 α 1 ) k α 0 ∈ B S . Note that α 1 is not modified by 0-removal.

Dimension 2

There are two different removal operations (0 and 1-removal) for dimension 2.

′ 0 = 1(α 0 α 1 )α 0 = 4 ∈ B S 3.2.1 0-Removal It consists in removing a 0-cell C =< α 1 , α 2 > (b). Let B S = Cα 0 -C (Cα 0 is defined as above)
. This operation can be applied only if the following precondition is satisfied:

∀b ′ ∈ C, b ′ α 1 α 2 = b ′ α 2 α 1 .
This constraint corresponds, in the general case, to the fact that the degree of the vertex is equal to 2 (2 edges are incident to the vertex). If this constraint is not satisfied, we do not know how to join the cells incident to C, and it is then impossible to define the removal in a simple way. [START_REF] Elter | Etude de structures combinatoires pour la représentation de complexes cellulaires[END_REF] proposes a generalization of this operation, but it is complex and cannot be used for an automatic process, in particular in automatic image processing.

The G-map resulting from 0-removal is obtained by redefining α 0 for the darts of B S as follows:

∀b ′ ∈ B S , b ′ α ′ 0 = b ′ (α 0 α 1 ) k α 0 where k is the smallest integer such that b ′ (α 0 α 1 ) k α 0 ∈ B S .
Note that this redefinition of α 0 is the same as for dimension 1 but concerns different darts, since it is a 0-cell within a 2D object: cf. Fig. 3 (intuitively, in the general case, this operation consists in applying twice the 0-removal defined for dimension 1).

1-Removal

It consists in removing a 1-cell C =< α 0 , α 2 > (b). This can be achieved without any precondition. Let B S = Cα 1 -C. The resulting G-map is obtained by redefining α 1 for the darts of B S as follows:

∀b ′ ∈ B S , b ′ α ′ 1 = b ′ (α 1 α 2 ) k α 1 , where k is the smallest integer such that b ′ (α 1 α 2 ) k α 1 ∈ B S .
Examples of 1removal are presented in Fig. 4, and 5. For this last example, k = 2 since the removed edge is incident twice to the same vertex.

Dimension n

The general definition of i -cell removal for an n-dimensional G-map is an obvious extension of the previous cases. Let C be an i -cell to remove; when

′ 1 = 1(α 1 α 2 )(α 1 α 2 )α 1 = 4 ∈ B S (since 1(α 1 α 2 )α 1 ∈ B S ,
this dart belongs either to C and to Cα 1 ).

i < n -1 the operation can be applied only when, informally 4 the degree of C is equal to 2 (a vertex incident to exactly two edges or an edge incident to two faces or a face incident to two volumes. . . ). The i -removal consists then in redefining α i for the darts of

B S = Cα i -C in the following way: b ′ α ′ i = b ′ (α i α i+1 ) k α i , where k is the smallest integer such that b ′ (α i α i+1 ) k α i ∈ B S .
We obtain so the general definition of the i -removal operation:

Definition 3 (i-cell removal) Let G = (B, α 0 , . . . , α n ) be an n-G-map, i ∈ {0, . . . , n-1} and C =<> N -{i} (b) be an i-cell, such that: ∀b ′ ∈ C, b ′ α i+1 α i+2 = b ′ α i+2 α i+1 . Let B S = Cα i -C, the set of darts i-sewn to C that do not be- long to C. The n-G-map resulting from the removal of this i-cell is G ′ = (B ′ , α ′ 0 , . . . , α ′ n ) defined by: • B ′ = B -C; • ∀j ∈ {0, . . . , n} -{i}, α ′ j = α j | B ′ ; 5 • ∀b ′ ∈ B ′ -B S , b ′ α ′ i = b ′ α i ; • ∀b ′ ∈ B S , b ′ α ′ i = b ′ (α i α i+1 ) k α i , where k is the smallest integer such that b ′ (α i α i+1 ) k α i ∈ B S .
Note that G ′ can contain only one n-cell, and may even be empty if G contains only one i-cell.

Contraction

Informally, i -contraction consists in contracting an i -cell into an (i -1)-cell. Contraction is the dual of the removal operation. Informally, the dual of a subdivision is a subdivision of the same space, in which an (ni)-cell is associated with each initial i -cell, and incidence relations are kept. A nice 4 The formal precondition is:

∀b ′ ∈ C, b ′ α i+1 α i+2 = b ′ α i+2 α i+1 . Note that if i = n-1
this condition doesn't apply and we can always remove any (n -1)-dimensional cell. 5 α ′ j is equal to α j restricted to B ′ , i.e. ∀b ∈ B ′ , bα ′ i = bα i property of G-maps is the fact that the dual G-map of G = (B, α 0 , . . . , α n ) is G ′ = (B, α n , . . . , α 0 ): we just need to reverse the involution order.

We can thus easily deduce the definition of i -contraction from the general definition of i -removal. We just have to replace '+' by '-' for indices of involutions for preconditions and operations, i.e. α i+1 α i+2 →α i-1 α i-2 and α i α i+1 →α i α i-1 .

(see two examples of contraction in Fig. 6 and7).

Definition 4 (i-cell contraction) Let G = (B, α 0 , . . . , α n ) be an n-G-map, i ∈ {1, . . . , n} and C =<> N -{i} (b) be an i-cell, such that 6 : ∀b ′ ∈ C, b ′ α i-1 α i-2 = b ′ α i-2 α i-1 . Let B S = Cα i -C , the set of darts i-sewn to C that do not belong to C. The n-G-map resulting from the contraction of this i-cell is G ′ = (B ′ , α ′ 0 , . . . , α ′ n ) defined by: • B ′ = B -C; • ∀j ∈ {0, . . . , n} -{i}, α ′ j = α j | B ′ ; • ∀b ′ ∈ B ′ -B S , b ′ α ′ i = b ′ α i ; • ∀b ′ ∈ B S , b ′ α ′ i = b ′ (α i α i-1 ) k α i , where k is the smallest integer such that b ′ (α i α i-1 ) k α i ∈ B S .

Generalisations

Previous definitions enable us to remove or to contract a single cell. For some applications, it could be more efficient to simultaneously apply several operations. Concretely, let G be an n-G-map. Assume that each dart belonging to a removed or contracted cell is marked with the dimension and type of the corresponding operation. Operations can be simultaneously applied if and only if:

• the cells are disjoined (involving that a dart has at most a unique mark and thus that there is no ambiguity when redefining α i for each i, 0 ≤ i ≤ n); • the preconditions of the corresponding operations are satisfied. Indeed, when a precondition is satisfied before a set of operations, this precondition is still satisfied after each step of the operation (the converse is obviously false), because cells are disjoined. This allows us to apply simultaneously or successively a set of operations and to obtain the same result.

We now present this generalization in several steps. First, we show that it is possible to simultaneously perform removals (resp. contractions) of several i -cells for a given i (0

≤ i ≤ n).
Generalisation 1 We can easily prove that the previous definition of removal (resp. contraction) stands for the removal (resp. contraction) of a set of cells of same dimension i. The (possible) precondition of the initial operation has to be satisfied for each cell (cf. Fig. 8(a)).

Moreover, removing (resp. contracting) simultaneously several i -cells or applying successively, in any order, the initial operation for any removed cell, produce the same result. The main idea of the proof is: each

α i redefinition con- sists in: ∀b ∈ B Si , bα ′ i = b(α i α i+1 ) k α i .
The darts of this path can be partitioned depending on the removed cells they belong to, i.e. bα

′ i = b(α i α i+1 ) k 1 (α i α i+1 ) k 2 . . . (α i α i+1 ) kp α i .
Each subpath corresponds to a single removal and so order is not important, and each removal does not depend on other operations.

We retrieve here the notion of connecting walk of Brun and Kropatsch [START_REF] Brun | Pyramids with combinatorial maps[END_REF]. Darts of B S i are surviving darts and α ′ i put in relation two darts of B S i by passing through some non-surviving darts (darts crossed by (α i α i+1 ) k ).

Generalisation 2

The previous generalization can be directly extended for simultaneously removing and contracting cells of same dimension i. A cell is either removed or contracted, but not both at the same time. The (possible) precondition of the corresponding initial operation has to be satisfied for each cell (cf. Fig. 8(b)).

More precisely, let CS i (resp. CC i ) be a set of i -cells to remove (resp. contract), such that CS i ∩CC i = ∅ and such that the (possible) precondition of i -removal (resp. i -contraction) operation is satisfied for each cell of CS i (resp. CC i ). Let satisfied for each cell. This last generalization is possible because the set of cells are disjoined. If we consider two disjoined cells, there are two possibilities: the two cells have the same dimension and they can be adjacent or not; the two cells have different dimensions and they can not be incident. The first case is covered by the previous generalizations. The second one is easy to consider since the cells are not incident: so some surviving darts exist between the two cells, and α i is redefined only for these darts. The other cells are not modified and so other preconditions are still valid when applying any subset of these operations.

B Si = (CC i ∪ CS i )α i -(CC i ∪ CS i ). As before, α i is redefined for these darts: ∀b ∈ B Si , bα ′ i = b ′ = b(α i α k 1 ) . . . (α i α kp )α i where p is the smallest integer such that b ′ ∈ B Si and ∀j, 1 ≤ j < p, if b c = b(α i α k 1 ) . . . (α i α k j-1 )α i ∈ CS i then k j = i + 1 else (b c ∈ CC i ) k j = i -1.

Generalisation 3 The previous generalization can be directly extended for the removal and/or contraction of a set of disjoined cells of any dimension. The (possible) precondition of the corresponding initial operation has to be

6 A pyramidal data structure for a spatial modeler

Introduction

For a few years we are interested in the design of a discrete modeling software [START_REF] Andres | Spamod: design of a spatial modeling tool[END_REF]. With this modeler, our goal is to create, import and handle different representations of a same geometrical object (mainly discrete and continuous representations). Modeling softwares exist that can handle discrete representations or Euclidean representations but rarely both. There are however situations where both representations could be useful.Let us give two examples to illustrate the long term goals of such a modeling system: Let us consider the general problem of computer vision as a first example. It is clear now that computer vision (for surveillance and security system applications for instance) requires low-end image processing techniques but also the construction of a high-level semantic model of the image content. There lies an important bottleneck in defining and constructing such a virtual high-level model. It would seem natural that a high-level model shouldn't be unrelated to virtual reality models as proposed in "continuous modeling". Being able to have both a continuous "high" level representation and a discrete "low" level one could be an asset.

As a second example, let us consider the problem of augmented reality with an image containing real world objects (acquired by photography for instance) and virtual, synthetic, ones. Having both coexist in a single image means mixing geometrically the objects, mixing artificial and natural light and shadow effects, etc. Augmented reality has many applications in movie and advertising, virtual surgery, architecture, etc. One bottleneck lies in the virtual models that haven't been conceived to handle low-end pixel or voxel information. Again, a close collaboration between both worlds should bring new insight in such problems.

In a first release of our modeler, the discrete and Euclidean forms of an object were coexisting, but were not directly linked. Then, in case of modification, each representation had to be entirely recomputed in order to keep coherent, the different representations. However, most operations are local ones; the effect of a local modification of a representation will be a local modification of the other representation. This requires the existence of connections between the different representations of an object.

In order to solve these problems, we propose the use of a pyramidal structure. Each level corresponds to a particular representation of the object: Euclidean, discrete analytical, borders of regions and pixel representations. Each level is linked with the level above and below it. This ensures the coherence between all the representations.

We present here our structure and some basic operations. For the moment, we deal only with 2D objects, because the discrete analytical reconstruction of a Euclidean object from a discrete object is still under investigation in 3D, while the 2D case has been fully solved [START_REF] Breton | Towards an invertible euclidean reconstruction of a discrete object[END_REF]. Note that we can extend the pyramidal structure and the discretization process of a Euclidean object in any dimension.

Our structure: a pyramid of G-maps

Our pyramidal structure consists of four inter-dependent levels (see Fig. 9): -Level 0: the discrete level This level corresponds to the discrete matrix representation which is composed of pixels. A label, such as grey scale or color, is associated with each pixel.

-Level 1: the discrete border level The discrete border level corresponds to the contours obtained for each uniform labelled 4-connected region. These contours are the interpixel ones [START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF][START_REF] Khalimsky | Boundaries in digital planes[END_REF] (note that these contours are a Cellular Complex as defined by [START_REF] Alexandroff | Diskrete raüme[END_REF]). We can see on Fig. 9 that these contours are made up of 0-cell and 1-cell series.

-Level 2: the discrete analytical level This level corresponds to the discrete analytical description of the region's contours. More precisely, each contour is described as a discrete analytical polygon.

-Level 3: the Euclidean level In this level, each region is described as a Euclidean polygon.

For each level, there are structure and shape information. More precisely, each level is composed of a generalized map which represents the topological information. Specific geometrical information (called embeddings) are associated with specific orbits of each level in order to define the shape of the object.

New levels could be added without major upheaval in the structure. Note that the levels of the pyramid are chosen in order to obtain progressive evolution between the generalized maps of the first and the last levels: each intermediary map can be obtained from the previous one by minor modifications.

Links between levels

Since our goal is to quickly reflect a modification of a level of the pyramid on the other levels, we set up bidirectional connections between the different levels of the pyramid, and more precisely between darts of the levels (surviving darts). Figure 10 shows some connections between darts of the pyramid (this is only a partial representation). The advantage of these connections is that all pyramid's embeddings become accessible for all levels. The information redundancy inside the pyramid is thus limited. For example, color embeddings are only associated with level 0 map, and each level can access these information through existing connections.

Fig. 10. Some connections between levels.

Level 0: the discrete level

The G-map of this level corresponds to a space subdivision into pixels (see Fig. 11(a)). Each pixel is represented by a square topological face associated with a label embedding. Moreover, integer coordinates are attached to each topological vertex.

Level 1: the discrete border level

The G-map of this level corresponds to a space subdivision into regions (see Fig. 11(b)), in which each topological edge corresponds to a 1-cell. This map is obtained by merging pixels of level 0 having the same label, by applying 1-removals (see Sec. 7). There are no embeddings at this level since level 0 embeddings (integer coordinates of 0-cells) can be accessed by using the links between levels of the pyramid.

Level 2: the discrete analytical level

The G-map of this level corresponds to a space subdivision into regions (see Fig. 11(c)), in which each topological edge corresponds to a discrete analytical segment. The discretization model used here is the discrete analytical model

Fig. 11. The four levels of the pyramid.

called Standard Model [START_REF] Andres | Modélisation analytique discrète d'objets géométriques[END_REF][START_REF] Andres | Discrete linear objects in dimension n: the standard model[END_REF]. We use this model because this is the only one that allows the analytical discretization of 4-connected polygons (segments and points). Moreover, this model is defined for any dimension. We can see in Fig. 12 an example of a standard line and its analytical description. A discrete standard line of parameters a, b, c is the set of integer points (x, y) verifying -ω ≤ ax + by + c < ω with ω = |a|+|b| 2 and a ≥ 0 or a = 0 and b ≥ 0. Level 2 map is obtained from level 1 map by applying 0-removals (see Sec. 7).

The embeddings of this level consist in a set of inequations. Each discrete segment is described by two opposite inequations associated with a topological edge, and each discrete vertex is described by four inequations associated with a topological vertex. To ensure the coherence between this level and the two previous ones, analytical inequations of level 2 must be satisfied by all 0-cells of level 1 (pixel centers of level 0).

Level 3: the Euclidean level

There are many different works [START_REF] Bertrand | Spécification algébrique et réalisation d'un modeleur interactif d'objets géométriques volumiques[END_REF][START_REF] Bertrand | Algebraic specification of a 3d-modeler based on hypermaps[END_REF][START_REF] Dufourd | Functional specification and prototyping with combinatorial maps[END_REF][START_REF] Vidil | Développement d'un modeleur à base topologique[END_REF] about modeling of Euclidean objects with maps. The G-map of this level corresponds to a generalized map with vertex embeddings that are 2D Euclidean coordinates (see Fig. 11(d)). The standard discretization of this level must correspond to the level 1.

Note that this level is not unique. Indeed, a discrete object can be seen as an equivalence class of Euclidean objects which have the same discretization. In our software we consider only one member of this class. It is thus not possible to ensure that, if we discretize a Euclidean object and then reconstruct it, we recover the same Euclidean object. On the other hand, a discrete object that is reconstructed, then discretized will remain unchanged. This supposes that the reconstruction process is deterministic. Since it depends on the starting point, a starting point choice convention is set up. Basically, the discrete point that has the smallest x value and among those the one with the smallest y value is the starting point.

Construction steps between consecutive levels

In this section, we present all algorithms we use to build a level from a neighbor level. Indeed, depending on the application, we want to build a level i given a level i+1 (discretization) or in the other way build a level i+1 given a level i (continuation).

Conversions between level 0 and level 1

• Level 1 is obtained by merging all incident faces of level 0 having same label. This corresponds to the use of 1-removal operations (see Sec. 3). The algorithm is the following: first, all edges between two pixels having the same label are marked. Then, these edges are removed during a single scan by applying simultaneous 1-removals (see Sec. 5). We can see on Fig. 13 an example of discrete border level, as well as some connections between the two levels. Note that each border dart of level 0 (i.e. without adjacent dart) is connected with one dart of level 1. Integer coordinates embeddings of level 0 can thus be accessed by level 1 through these connections. • To build level 0 from level 1, we need to reconstruct the pixels, for example by using a flood-fill algorithm on a matrix of 1-cells corresponding to level 1. We can see in Fig. 14(b) the level 0 obtained from Fig. 14(a). Note that the discrete space represented in Fig. 14

Conversions between level 1 and level 2

• To obtain level 2 from level 1, we use an operation called analytical recognition [START_REF] Breton | Towards an invertible euclidean reconstruction of a discrete object[END_REF]. During the analytical recognition step, all vertices that belong to the same discrete segment are removed (by using 0-removal operation) in order to keep only the two edge extremities. We can see in Fig. 15 the simplification made after a discrete segment recognition. At the end of the process, we have recognized all edges of the initial object. We can see in Fig. 16 the final result obtained for level 2. The generalized map has been simplified and each topological edge corresponds to a recognized discrete segment. Analytical inequations computed during the analytical recognition are linked with each edge and with each vertex of the map. • To build level 1 starting from level 2, we need to discretize the contours of level 2 (the discretization model used is the standard model, see Sec. 6). For each edge of level 2, if the standard discretization of the two edge end points differ from each other (see Fig. 17), the edge is incrementally discretized from one end point to the other, and a new vertex is inserted into the edge while the other end point is not reached. If the standard discretization of the two edge end points are equal (see Fig. 18), all darts of the edge are marked. Then, all marked edges are contracted using simultaneous 1-contractions (see Sec. 5).

However, some problems can arise such as degenerated faces (see Fig. 19). These faces disappear in the discrete level due to the loss of information in- herent to the discretization process. Thus, we chose to delete these faces mainly using 2-contractions (see Sec. 4).

Conversions between level 2 and level 3

Level 2 is topologically quite similar to level 3, but the main difference between these two levels consists in embeddings. Indeed, level 2 has analytical inequation embeddings whereas level 3 has Euclidean coordinate embeddings.

• To build level 3 from level 2, we make a copy of level 2, and compute the Euclidean coordinate embeddings during the analytical reconstruction step (see [START_REF] Breton | Towards an invertible euclidean reconstruction of a discrete object[END_REF] for more details on particular cases that can occur during the analytical reconstruction).

• When we build level 2 starting from level 3, we also start with a copy of level 3, but we can eventually make some simplifications. For example, the gray square in Fig. 20(a) shows a set of vertices having the same standard discretization. All these vertices, except the two extremities, can thus be removed using 0-removals. We keep the two extremities in order to keep the inequations of the other edges, that are computed from the Euclidean coordinates.

Conclusion

In this paper, we have defined removal and contraction operations, that can be applied to any cells of any G-maps, whatever their respective dimensions. Moreover, we have studied how to simultaneously perform different operations. These definitions are homogeneous for any dimension. Since combinatorial maps [START_REF] Braquelaire | A topological structuring for aggregates of 3d discrete objects[END_REF][START_REF] Brun | Dual contraction of combinatorial maps[END_REF] can be easily deduced from orientable generalized maps [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasimanifolds[END_REF], these operations can also be defined on combinatorial maps.

We intend to revisit the works of Brun and Kropatsch for handling irregular pyramids of n-dimensional generalized maps. Properties of removal and contraction operations would enable us to establish relations between two contiguous levels within a pyramid, and thus between any levels. Efficient data structures could be deduced, taking these relations into account.

We have shown in a first experiment that we can easily define a combinatorial pyramid, thanks to the basic operations. This is a 2D pyramid, but it can be extended for higher dimensions, due to the general definition of removal and contraction operations for any dimension.

Moreover, this experiment shows a new way for using combinatorial pyramids. Indeed, a pyramid is generally used in order to represent an object at different levels of details. Here, we handle different representations of a same object, from continuous to discrete representation.This allows us to work with the best representation for a given operation. Moreover, we can progressively modify our structure, through the different intermediate levels. This simplifies the algorithms and the operations which maintain the coherence between the different levels of the pyramid. Now we are studying the definition of several operations, such as modification or consultation algorithms. Classical Euclidean operations can be translated into discrete operations, by discretizing level 3. We are now considering typical discrete operations that will be translated into Euclidean operations by analytical reconstruction and analytical recognition. The main goal of our work is now to study the way of updating a pyramid after a given operation, in order to maintain the coherence between the different levels.

Fig. 1 .

 1 Fig. 1. (a) A 2D subdivision. (b) The corresponding 2-G-map (involutions are given explicitly in the array). Darts are represented by numbered black segments. Two darts in relation by α 0 share a little vertical segment (ex. darts 1 and 2). Two darts in relation by α 1 share a same point (ex. darts 2 and 3). Two distinct darts in relation by α 2 are parallel and close to each other (ex. darts 3 and 20); otherwise, the dart is its own image by α 2 (ex. dart 2). Dart 1 corresponds to (s 1 , a 1 , f 1 ), dart 2 = 1α 0 corresponds to (s 2 , a 1 , f 1 ), 3 = 2α 1 corresponds to (s 2 , a 2 , f 1 ), and 20 = 3α 2 corresponds to (s 2 , a 2 , f 2 ). The vertex incident to dart 2 is < α 1 , α 2 > (2) = {2, 3, 20, 21}, the edge incident to dart 3 is < α 0 , α 2 > (3) = {3, 4, 19, 20}, and the face incident to dart 9 is < α 0 , α 1 > (9) = {9, 10, 13, 14, 17, 18}.

Fig. 2 .

 2 Fig. 2. 0-removal in 1D. (a) Initial 1-G-map. (b) Result. C =< α 1 > (2) = {2, 3} (darts marked with empty squares), Cα 0 = {1, 4} = B S (darts marked with crosses). 0-removal consists in setting 1α ′ 0 = 1(α 0 α 1 )α 0 = 4 ∈ B S and 4α ′ 0 = 4(α 0 α 1 )α 0 = 1 ∈ B S .

Fig. 3 .

 3 Fig. 3. 0-removal in 2D. (a) Initial 2-G-map. (b) Result. C =< α 1 , α 2 > (2) (darts marked with empty squares), Cα 0 = B S (darts marked with crosses). For instance, 1α ′ 0 = 1(α 0 α 1 )α 0 = 4 ∈ B S

Fig. 4 . 1 -

 41 Fig. 4. 1-removal in 2D in the general case. Darts of the edge to remove are marked with circles. (a) Initial 2-G-map. (b) Result.

Fig. 5 . 1 -

 51 Fig. 5. 1-removal in 2D of a loop. (a) Initial 2-G-map. (b) Result. For instance, 1α′ 1 = 1(α 1 α 2 )(α 1 α 2 )α 1 = 4 ∈ B S (since 1(α 1 α 2 )α 1 ∈ B S ,this dart belongs either to C and to Cα 1 ).

Fig. 6 . 1 -

 61 Fig. 6. 1-contraction in 1D. Darts of the edge to contract are marked with black disks. (a) Initial 2-G-map. (b) Result.

Fig. 7 . 1 -

 71 Fig. 7. 1-contraction in 2D. (a) Initial 2-G-map. (b) Result.

Fig. 8 .

 8 Fig. 8. (a) Simultaneous 1-removals in 2D. (b) Simultaneous 1-removals (dart marked with circles) and 1-contractions (dart marked with black disks) in 2D. For instance, 1α ′ 1 = 1(α 1 α 0 )(α 1 α 2 )α 1 = 2 since the edge incident to 1α 1 is contracted and the edge incident to 1(α 1 α 0 )α 1 is removed.

Fig. 9 .

 9 Fig. 9. The pyramidal structure.

5 Fig. 12 .

 512 Fig. 12. Example of standard line and its analytical description. Pixels in gray correspond to the standard discretization of (b), i.e. they are pixels crossed by the line. (a) and (c) give the two opposite inequations of the analytical description. Note that the inequation (c) is a strict one. These inequations are satisfied by the pixel centers.

Fig. 13 .

 13 Fig.[START_REF] Breton | Towards an invertible euclidean reconstruction of a discrete object[END_REF]. Example of discrete (level 0) and border discrete (level 1) levels. (a) Level 0 map with marked darts (marks are represented by black circles). (b) Level 1 map obtained after several 1-removals. Arrows show some connections between darts, and all black darts are connected ones.

  (a) and Fig. 14(b) are not the same, but are duals.

Fig. 14 .

 14 Fig. 14. (a) A discrete border level (level 1) and (b) its associated discrete view (level 0).

Fig. 15 .

 15 Fig. 15. Building level 2 from level 1. (a) Black darts belong to the recognized discrete segment. (b) Marked darts are removed using 0-removals. (c) Result.

Fig. 16 .

 16 Fig. 16. Example of discrete border (level 1) and discrete analytical (level 2) levels. (a) Level 1 with marked darts. (b) Level 2 obtained after several 0-removals. Upper arrows show some connections between darts of level 1 and level 2.

Fig. 17 .Fig. 18 .

 1718 Fig. 17. First case. (a) The black edge is the processed one. Each filled pixel is crossed by the corresponding Euclidean segment. (b) Incremental discretization of the edge.

Fig. 19 .

 19 Fig. 19. Example of degenerated faces. (a) Euclidean object. (b) Map corresponding to the object in (a). (c) Standard discretization of the two black edges in (a). Black darts are contracted using two 2-contractions. (d) Standard discretization of the gray filled face in (a). Black darts are contracted using one 2-contraction.

Fig. 20 .

 20 Fig. 20. Simplification of the discrete analytical level. (a) Euclidean object. (b) The darts marked with black squares are removed by 0-removals. (c) Discrete analytical level. Note that all black darts are connected ones.

  ). Darts are represented by numbered black segments. Two darts in relation by α 0 share a little vertical segment (ex. darts 1 and 2). Two darts in relation by α 1 share a same point (ex. darts 2 and 3). Two distinct darts in relation by α 2 are parallel and close to each other (ex. darts 3 and 20); otherwise, the dart is its own image by α 2 (ex. dart 2). Dart 1 corresponds to (s 1 , a 1 , f 1 ), dart 2 = 1α 0 corresponds to (s 2 , a 1 , f 1 ), 3 = 2α 1 corresponds to (s 2 , a 2 , f 1 ), and 20 = 3α 2 corresponds to (s 2 , a 2 , f 2 ). The vertex incident to dart 2 is < α 1 , α 2 > (2) = {2, 3, 20, 21}, the edge incident to dart 3 is < α 0 , α 2 > (3) = {3, 4, 19, 20}, and the face incident to dart 9 is < α 0 , α 1 > (9) = {9, 10, 13, 14, 17, 18}.

The discretization model chosen here is the discrete analytical model called Standard Model.

An involution f on S is a one to one mapping from S onto S such that f = f -1 .

Let {Π 0 , . . . , Π n } be a set of permutations on B. The orbit of an element b relatively to this set of permutations is < Π 0 , . . . , Π n > (b) = {Φ(b), Φ ∈< Π 0 , . . . , Π n >}, where < Π 0 , . . . , Π n > denotes the group of permutations generated by Π 0 , . . . , Π n .

Note that this condition doesn't apply for i = 1, so we can always contract any edge.


Konrad Hinsen
email: kon-rad.hinsen@cnrs.fr.

Reusable vs. re-editable code

DEPARTMENT: Scientific Programming

Reusable vs. re-editable code

Konrad Hinsen Centre de Biophysique Moléculaire, Orléans One of my current projects is a MOOC (Massively Open Online Course) on reproducible research that I am working on within a multidisciplinary team. If all goes well it will go on-line in autumn. French-speaking readers might want to watch out for it on the FUN platform (https://www.fun-mooc.fr/). But my topic for today is not the MOOC, it's a lesson learned during its preparation. To sum it up in a keyword, it's about the good and bad sides of reusable code.

One module of the MOOC shows how to prepare a simple data analysis as a reproducible computational document. We decided to present this example in several variants, using different languages and technologies. One variant uses a Jupyter notebook and the Python 3 language. The second one is done in R, using RStudio as an authoring tool for the R Markdown format. The final variant, for power users, mixes several programming languages using the org-mode extension to the Emacs editor. Our goal is to demonstrate the current best practices for each platform. Rather than developing a single strategy for the data analysis and translating it into three dialects, we start from scratch three times, asking ourselves how a competent user of each platform would approach the problem. An additional benefit of this approach is the possibility to compare the strategies that the three platforms encourage, and that is what I will do in the following.

Since we concentrate on reproducibility rather than on mathematics or statistics, the chosen analysis is simple: find the distribution of the severity of the annual flu epidemics in France from data on the number of diagnosed cases for each week since 1985. The data is publicly available from the "Réseau Sentinelles" (http://www.sentiweb.fr/). The essential steps are (1) download the data as a CSV file, (2) read the CSV file, (3) check for potential problems, in particular missing data points, (4) transform the weekly sampling into an annual sampling, and (5) produce a histogram of the annual number of flu victims.

Python and pandas

Since I am the Python guy of the team, my first assignment was to prepare the analysis as a Jupyter notebook using the Python 3 language. The standard tool in the Python ecosystem for working with tabular data is the pandas library (http://pandas.pydata. I have shortened the column name line a to make it fit into the page. The meaning of the columns is documented on the "Réseau Sentinelles" Web site. We only need two of them: week and inc. The former encodes a reference to each week in a somewhat uncommon format: it's a six-digit string with the first four digits indicating the year and the last two the week number according to the ISO-8601 standard. The inc column contains the number of reported flu cases during each week.

Reading such a dataset using pandas isn't difficult:

import pandas as pd def data_url(year1, week1, year2, week2): return "http://websenti.u707.jussieu.fr/sentiweb/" + \ "api/data/rest/getIncidenceFlat?indicator=3" + \ ("&wstart=%4d%02d&wend=%4d%02d" % (year1, week1, year2, week2)) + \ "&geo=PAY1&$format=csv" data = pd.read_csv(data_url(1990, 1, 1995, 1), skiprows=1)

The argument skiprows=1 tells pandas to ignore the first line of the downloaded file, which is a comment that doesn't fit into the CSV format. What we get back as data is a DataFrame object with 262 rows and 10 columns. When reading a CSV file, pandas tries to guess appropriate data types for each column. In this case, it decides on "int64" for both of the columns we need, although this isn't quite appropriate for week. An excellent opportunity to explain the subtleties of data representations and their conversions to the students in our MOOC! Unfortunately, there is another subtlety to come. After testing the analysis code on the rather small five-year sample, we want to run it on the full dataset, which starts in 1985 and is continually updated. Using data url(1985, 1, 2018, 1), we get data for 33 years. On the screen it looks much the same. But a more careful inspection shows that the data type for inc has changed from "int64" to "object". Before going through the 1723 lines to check for an anomaly that could explain this change, I took a look at the pandas documentation for CSV import (http://pandas.pydata.org/pandasdocs/stable/io.html). It mentions "intelligent conversion of tabular data" but the rules for this intelligent action are not given. In view of the length of the documentation, which spreads over several pages, I concluded that going through the 1723 lines of data is the lesser evil. And indeed the offending line clearly stands out visually:

198919,3,-,,,-,,,FR,France A missing data point for week 19 of the year 1989 is indicated by a dash instead of a number in the inc column. Faced with a column containing many integers and one string, pandas apparently decided on the "object" type.

I won't discuss the many options of dealing with such nuisances here. They are the rule rather than the exception when dealing with data coming from the real world. In fact, the pandas documentation I just cited proposes a few optional arguments I could pass to the function read csv in order to deal with the irregularity in the input file. What makes this problem unpleasant is not that it's difficult to deal with, but that it's difficult to identify. A data type chosen as a function of the data values means that the subsequent analysis code can never be trusted to work in general. All you can say is that it works for the data ranges that it has been explicitly tested on.

Emacs and org-mode

I won't say much about the R version of the data analysis, since it was prepared by my R-savvy colleagues. It isn't very different in spirit from the Python 3/pandas version, and requires similar decisions for how to deal with irregularities. What I will focus on is the multi-language version in the Emacs org-mode environment, which has been presented in the Scientific Programming department before [START_REF] Schulte | Active Documents with Org-Mode[END_REF] . In org-mode, so-called "source code blocks" contain executable pieces of code, much like the cells in a Jupyter notebook. However, org-mode offers much more flexibility in combining these source code blocks into a computational document. Their execution can be in sequential order, as for a Jupyter notebook, but also in arbitrary order, defined by naming the inputs and outputs of each block and thereby defining a a data-flow graph. But the main difference for this story is that each source code block can be written in a different language, which makes exchanging data between them much more challenging. Org-mode defines a common data representation that it translates to and from for each language. This common representation is a tree structure made from nested lists, with the leaves being numbers or strings. Such tree structures are well-known to be very flexible, and therefore widely used in data interchange formats, the best-known one probably being XML.

The requirement to use this tree data structure at the interface between source code blocks discourages the use of language-specific data structures such as pandas' DataFrames, which would require a conversion layer at the beginning and end of each snippet of code. What is encouraged instead is to keep all data in the tree structure as much as possible. This corresponds to a long tradition of data representation in the Lisp family of programming languages, which is no coincidence because most of Emacs, and in particular all of org-mode, is written in Emacs' own Lisp dialect called Emacs Lisp. This emphasis on a single universal data structure, in stark contrast with the ideas of object-oriented programming, is often justified using a quote from computer science pioneer Alan J. Perlis [START_REF] Perlis | Epigrams on Programming[END_REF] dating back to 1982: "It is better to have 100 functions operate on one data structure than 10 functions on 10 data structures."

Emacs being a text editor, it has excellent support for manipulating text, and therefore it seems natural to start the flu data analysis using Emacs Lisp. A side benefit is the possibility of loading the data into an Emacs buffer, which allows direct inspection using all of Emacs' text editing functions. The first step thus is to download the data into a read-only buffer named * data * :

(require 'url) (defun data-url (year1 week1 year2 week2) (concat "http://websenti.u707.jussieu.fr/sentiweb/" "api/data/rest/getIncidenceFlat?indicator=3" (format "&wstart=%4d%02d&wend=%4d%02d" year1 week1 year2 week2) "&geo=PAY1&$format=csv")) (with-current-buffer (get-buffer-create " * data * ") (unless buffer-read-only (url-handler-mode) (insert-file-contents (data-url 1990 1 1995 1)) (setq buffer-read-only t)))

Making the buffer read-only provides a protection against accidental manipulation. Since the data format is actually very simple, I wrote my own special-purpose parser for this particular file rather than searching for a CSV parser library for Emacs Lisp. All it takes is (1) split the file into lines, (2) discard the first line, (3) split each of the remaining lines at the commas, yielding the columns of the table, and (4) select the first and third columns. In Emacs Lisp this can be written as (require 'dash) (with-current-buffer " * data * " (let * ((lines (split-string (buffer-string) "\n" t)) (data-lines (rest lines)) (table (--map (split-string it ",") data-lines))) (-insert-at 1 'hline (-select-columns '(0 2) table ))))

The last operation, inserting a hline between the column names and the data, puts the parsed data into the format of an org-mode table, which Emacs then formats nicely before inserting it into the computational document. For the benefit of readers not familiar with Emacs Lisp, and for a better comparison with the pandas-based code I have shown earlier, here is an implementation of the same method in Python 3. It leads to exactly the same result, but doesn't permit inspection of the data file in an Emacs buffer. Note that the final return is a particularity of Python embedding in org-mode: a Python source code block is treated as the body of a function. The None being inserted after the column header is the Python equivalent of the hline in Emacs Lisp.

Not counting the shared function data url, this back-to-the-basics implementation consists of 10 lines of code, compared to two lines for reading the data using pandas. However, and that is the main message of this story, it is immediately clear to anyone with a minimal knowledge of Python 3 how the code works, which assumptions it makes about the format of the data, and that all the data is always returned as text strings. There is no hidden intelligence in these 10 lines of basic Python 3 code.

Another way to compare the two Python 3 implementations is to compare the effort that a reader must invest to understand what is going on. The pandas version provides an immediate "executive summary": we are reading a CSV file. Moving on to a more detailed understanding of what we get in return, and under which conditions we can expect the code to work, is almost a research project, since even the multi-page documentation of the read csv function does not provide the answers. The basic Python 3 version, on the other hand, requires more effort to reach a first level of understanding, because the only superficially available information is that the code processes a text file. However, after a two-minute inspection of the code, a reader knows everything there is to know, without reading lengthy documentation or doing experiments.

As I was working on this example for the MOOC, I stumbled upon a ten-year-old interview with Donald Knuth (http://www.informit.com/articles/article.aspx?p=1193856) that I saw quoted in the context of literate programming, one of Knuth's major inventions. Since time spent reading Donald Knuth is rarely wasted, I read through the complete interview and found the following statement: I also must confess to a strong bias against the fashion for reusable code. To me, "re-editable code" is much, much better than an untouchable black box or toolkit. I could go on and on about this. If you're totally convinced that reusable code is wonderful, I probably won't be able to sway you anyway, but you'll never convince me that reusable code isn't mostly a menace. This is a nice summary of my experience preparing two versions of my data analysis example. The read csv function of pandas is designed to be reusable in a wide range of situations, but that also makes it a nearly impenetrable black box. In contrast, the much more basic approach that the org-mode conventions encouraged me to adopt is not reusable but re-editable. People who need to process a somewhat different dialect of the CSV format can read it, understand it, and adapt it to their own situation.

And the winner is...

It may seem at this point that I am fully adhering to Donald Knuth's preference for reeditable over reusable code, but it is not as simple as that. If I have explained the advantages of re-editable code in so much detail, it is mainly because I have rarely seen them exposed. It is much easier to find discussions of the virtues of reusable code. The two main reasons given are economy and reliability. The economy argument is based on development costs being shared among many users. This does, however, silently suppose that the development investment leads to savings for the users, which isn't obvious as my CSV reading example illustrates. The reliability argument says that heavily used code is better tested, an idea often referred to as "Linus's law" in the open-source developer community and summarized as "given enough eyeballs, all bugs are shallow". It sounds plausible but, like most statements about software engineering, is hard to back up by solid evidence. Heavy usage has certainly not prevented the existence of a severe security bug in OpenSSL for two years [START_REF] Durumeric | The Matter of Heartbleed[END_REF] .

What, then, is the better approach: reusable or re-editable code? As so often, it depends. Let's look at the other end of my reproducible data analysis example: plotting a histogram of the annual flu incidence. I wouldn't even think of doing this in basic Python, or whatever other language, because of the complexity of the task. Even the people who write reusable plotting libraries don't limit themselves to basic Python or C, but build on reusable graphics libraries. So what's the difference between reading a CSV file and plotting a histogram? It's the localization of complexity.

In the CSV reading case, all the complexity is at the surface. In fact, it's the CSV format itself that is the source of complexity. Or perhaps one should better say that it is the habit of referring to a diverse family of file formats by a single name, leading to the idea that a single piece of code should deal with all of them. This complexity cannot be hidden from the user, which is why the documentation of pandas' read csv function is so long. In the histogram plotting case, most of the complexity is in the depth of the implementation. The user of a plotting library does not need to know how graphics libraries work or how the PDF format is defined, and would in general prefer not to be exposed to such details.

As a rough guide to whether packaging a piece of code as a reusable library makes sense, or whether looking for reusable code is preferable to rolling your own code, compare the length of the code to the length of documentation that users would consider sufficient. Reusable code makes sense when the former clearly exceeds the latter. Another way to look at this is that well-written code is its own best documentation. Hiding the code from its users and replacing it by documentation makes sense only when this clearly reduces the users' mental load.

A corollary of this principle is that if you publish your code, make it either re-editable, i.e. self-explaining, or reusable, i.e. well documented. Scientists often complain that writing documentation is a lot of work that is not well rewarded. Going for re-editable and keeping documentation to a minimum is thus worth considering as an alternative, keeping in mind that the audience for reusable and re-editable code is not quite the same. And that is perhaps the ultimate criterion for choosing the right approach: the one that best suits your audience.

In the context of scientific programming, there is one more criterion to keep in mind. If your code implements a scientific method, rather than a technical detail such as reading CSV files or plotting histograms, the decision to expose or to hide the code implies the decision to expose or hide the method itself. Well-known and frequently used methods are the perfect candidates for implementation in reusable and documented libraries. Recent methods, and in particular methods still under development, are better proposed as reeditable implementations. Otherwise, readers might feel cheated when presented with a computational solution to a scientific problem that consists only of a call to a black-box library -for an illustrative example, see the public review of a submission to the ReScience journal at https://github.com/ReScience/ReScience-submission/pull/46#issuecomment-366257844. The solution in that case was to include the library code into the review, which also shows that reusable and re-editable code are just the extremes of a continuum.

  from urllib.request import urlopen handle = urlopen(data_url(1990, 1, 1991, 1)) data_as_bytes = handle.read() data = data_as_bytes.decode('ascii') lines = data.strip().split('\n') data_lines = lines[1:] table = [line.split(',') for line in data_lines] selected_columns = [(row[0], row[2]) for row in table] selected_columns.insert(1, None) return selected_columns

  org/), so that's what I based my code on. It starts by downloading and parsing the data, which looks like this: # reseau Sentinelles, INSERM, UPMC, http://www.sentiweb.fr week,indicator,inc,inc_low,inc_up,inc100,inc100_low, ... 199501,3,26556,21916,31196,46,38,54,FR,France 199452,3,22036,17495,26577,38,30,46,FR,France 199451,3,26912,21750,32074,47,38,56,FR,France 199450,3,28831,23470,34192,50,41,59,FR,France ...