T Caissard

D Coeurjolly

J O Lachaud

T Roussillon

Laplace-Beltrami operator on Digital Curves

Keywords: JFIG 2016, laplacians, heat equation, convolution, digital surfaces

Many problems in image analysis, digital processing and shape optimization are expressed as variational problems and involve the discritization of laplacians. Indeed, PDEs containing Laplace-Beltrami operator arise in surface fairing, mesh smoothing, mesh parametrization, remeshing, feature extraction, shape matching, etc. The discretization of the laplace-Beltrami operator has been widely studied, but essentially in the plane or on triangulated meshes. In this paper, we propose a digital Laplace-Beltrami operator, which is based on the heat equation described by [BSW08] and adapted to 2D digital curves. We give elements for proving its theoretical convergence and present an experimental evaluation that confirms its convergence property.

Introduction

Many problems in image analysis, digital processing and shape optimization are expressed as variational problems and involves the discritization of laplacians. Indeed, PDEs containing Laplace-Beltrami operator arise in surface fairing, mesh smoothing, mesh parametrization, remeshing, feature extraction, shape matching, etc (see for example [START_REF] Lévy | Spectral Mesh Processing[END_REF]). For instance, computing geodesics on triangulated surfaces can be formulated from the heat diffusion equation [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF] or by a limit probability distribution of Dirac measures in optimal transportation problems [START_REF] Solomon | Earth Mover's Distances on Discrete Surfaces[END_REF].

The subject of the laplace-Beltrami operator and its discretization has been widely studied. Some approaches rely on the theory of exterior calculus which has been developped in the computational mathematics and geometry processing community, with focus on triangular meshes. Purely combinatory approches can be found in [START_REF] Grady | Discrete calculus: Applied analysis on graphs for computational science[END_REF]. Another approach developped in [START_REF] Hirani | Discrete exterior calculus[END_REF] relies more on the relation between the discrete structure and the continuous one. On triangular meshes, those formulations are closely related to the famous cotangent formula [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF], and is equivalent to specific cases of finite element method [START_REF] Lévy | Spectral Mesh Processing[END_REF].

An important objective when proposing discretizations of the laplace-Beltrami operator is to give convergence results: as meshes refine and tend toward the underlying manifold under certain properties, the approximated laplace-Beltrami operator should tends toward the usual one on the manifold. On arbitrary triangular meshes, it is shown that the computed laplace-Beltrami operator cannot recover all the properties of the smooth manifold one [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF]. Regard-ing discrete exterior calculus, Hildebrandlt et al. [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF] provided convergence results when the triangulated meshes tend toward the continuous manifold with those properties: Hausdorff distance tends to zero, mesh normals tend to surface normals and the mesh is projected one-to-one on surface. Similar proofs exist in the context of finite element method [AFW06a, AFW06b] and for chainlet discrete calculus [START_REF] Harrison | Geometric hodge star operator with applications to the theorems of gauss and green[END_REF].

In the case of triangular meshes, as the structure becomes thinner, estimated intrinsic geometric informations converge toward the real ones and can be easily injected inside the calculus. In our case, we want to define the Laplace-Beltrami operator on a different geometric structure: it is called a digital surface and is the border of a subset of the digital space Z d . One can see them geometrically as unions of faces of d-dimensional cubes. Digital surfaces can be defined as digitizations of manifolds but such digitizations present several difficulties: points are spaced evenly, they do not interpolate the smooth surface and elementary normals are not informative. Therefore, in opposition to triangular meshes, one needs to be more clever to recover geometric informations such as normals, areas or curvature for example. Fortunately, this topic was studied in depth in the last ten years and several progresses were made during this period. Parameter-free tangent and normal estimation along 2D digital curves were established [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF][START_REF] De | Convex digital polygons, maximal digital straight segments and convergence of discrete geometric estimators[END_REF] using properties of maximal digital straight segments. This approach was extended to 3D curves in [START_REF] Postolski | Tangent estimation along 3d digital curves[END_REF]. Further works using digital version of integral invariants [CLL14, PWHY09] induced convergent estimation of the normal vector field along digital surfaces in arbitrary dimension as well as the whole curvature tensor [CLL14, LCL14]. As shown on Fig. 1 and Fig. 5, using classical DEC operator, one cannot construct a convergent Laplace-Beltrami operator in 1D even by using convergent geometric estimators. Indeed, theoretical analysis on 2D curves shows that one cannot expect convergence of the dicretized Laplacian using exterior calculus toward the real Laplace-Beltrami on Riemannian Manifold. Recent work [START_REF] Thibert | Properties of gauss digitized sets and digital surface integration[END_REF] studied the closeness between a digital surface in a space of dimension d and its underlying smooth manifold structure. They particulary proved convergence of discretized integrals using convergent estimated normals on digital surface.

In this paper, we address the problem of computing a digital Laplace-Beltrami operator en 2D digital curves by using the work of [BSW08] based on the heat equation. We focus our study on 1D curves embedded in 2D as even in this simple case, the discritization of the operator is not trivial. We give elements for proving its theoretical convergence and present an experimental evaluation that confirms its convergence property.

Digital Surfaces

In all the paper, we assume that there is some ideal shape M in the space with a smooth topological boundary. We recall the definition of the Gauss Digitization process, which makes the link between the continuous shape and its digital approximation:

Definition 1 Let h > 0 be the sampling grid step. The Gauss Digitization of a shape M ⊂ R d is defined as D h (M) := M ∩ (hZ) d , where d is the dimension of the space.

The digitization process has therefore a very simple scheme: it considers the discrete points of the infinite regular grid with sample rate h and keeps only the ones inside the shape (see Fig. 2). We need also the notion of h-cube of a discrete point: for some z ∈ (hZ) d it is the closed d-dimensional axisaligned cube of R d centered on z with edge length h. We denote it by Q h z . The collection of all h-cubes of the discrete points of a digitized shape is called its h-cube embedding. Finally, we arrive at the structure on which we will do computations:

Definition 2 The digitized h-boundary of M, denoted by ∂ h M, is the topological boundary of the h-cube embedding of the Gauss digitization of M:

∂ h M := ∂ z∈D h (M) Q h z .
Geometrically it is a "staircased" surface that approximates the boundary of the continuous shape M, and is also called a digital surface.

We wish to prove that our digital laplacian on the digitized boundary converge to the standard laplacian on the continuous shape boundary. The convergence setup will be the following: given an estimated operator on the digital surface and a real operator on the smooth structure, we look at pointwise convergence between those two quantities as the grid step h tends toward zero (that is when ∂ h M → h→0 ∂M in the Hausdorff sense). More details about digital surfaces and their relationship with the underlying manifold can be found in [START_REF] Thibert | Properties of gauss digitized sets and digital surface integration[END_REF].

Figure 2: Illustration of the notations used in this paper. A smooth shape M, its Gauss digitization D h (M), its h-cube embedding Q h D h (M) and its h-boundary ∂ h M

Discretization on triangular meshes

We first describe the work of [BSW08], which proposes a discrete laplacian for triangulated meshes. Let g 0 : ∂M → R be a smooth function (ie the initial temperature distribution on ∂M). Let g : ∂M × (0, T ) → R a function which solves the partial differential equation called the Heat Equation:

∆g = ∂ ∂t g, (1)
with inital condition g(0) := g(•, 0) = g 0 : ∂M → R. The function g(t) := g(•,t) is therefore the temperature distribution on M at time t. In R d , one can find (see [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF] for example) an exact solution:

g(x,t) = 1 (4πt) d 2 y∈R d e -||x-y|| 2 4t g 0 (y)dy, (2)
where

H t R d := 1 (4πt) d 2 e -||x-y|| 2 4t
is called the heat kernel. On an arbitrary manifold, the heat kernel involves a more complex construction (see [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF]). Fortunately, we know from [START_REF] Molchanov | Diffusion processes and riemannian geometry[END_REF] that on a compact subset of ∂M the real heat kernel is not too far from the R d one.

The very idea of [BSW08] is to make the heat equation not a property of the laplacian but a definition by injecting Eq.2 into Eq.1:

∆g(x, 0) = lim t→0 1 t y∈∂M H t ∂M (x, y)g(y)dy -f (x) (3a) = lim t→0 1 t y∈∂M H t ∂M (x, y)(g(y) -g(x))dy. (3b)
Then they give the following discretization of the laplacian at a vertex w of a triangular mesh K:

L t K g(w) := 1 4πt 2 ∑ p∈V e -||p-w|| 2 4t g(p) -g(w) A(p), ( 4
)
where V is the set of vertices of K and A represents the mean area of triangles that surround p. They prove a convergence theorem under the assumption that vertices of K lie on the underlying manifold, and some additional properties on K (see [BSW08]).

On digital structures, we must show how to compute both || • || and A. The former relies on computing convergent distances in 2D, and areas in 3D which is straightforward. The later becomes non-trivial in as one must know how to compute distance between two points on a surface. Therefore, as stated in the introduction, we present the discritization of the operator on 1D curves embedded in 2D.

Discretization on digital curves

We describe how to discritize || • || and A from Eq.4. First, we adapt Eq.4 onto our digital curve:

∆h g( ẇ) := 1 √ 4πt ∑ ċi∈∂hM e -d( ċi , ẇ) 4t g( ẇ)-g( ċi ) µ h ( ċi ), (5)
where ẇ and ċ are vertices (element of dimension 0 of ∂ h M as in Fig. 4) on our digital curve, d is the approximated distance between two vertices and µ h is the elementary measure associated to a vertex. To measure such lengths we use the Integral Invariant normal estimator [START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF], whose convergence speed is in O(h 2 3 ). For each linel l of the digitized boundary (a linel is an edge joining two vertices), we define its measure µ h (l) as:

µ h (l) := n(l), ne(l) , (6)
where •, • is the dot product, n the estimated normal and ne the trivial normal to the linel (see Fig. 3). It is known (see [START_REF] Coeurjolly | Multigrid Convergence of Discrete Geometric Estimators[END_REF]) that this measure gives a good approximation of the real projected arclength on the curve. Then we have:

µ h ( ċi ) = µ h (l i ) + µ h (l i+1 ) 2 and d( ċi , ẇ) = ∑ j µ h (l j ), (7)
where j lies in the range of indices corresponding to the linels between ċi and ẇ. The µ h function is the digital analog of A in Eq.4 on triangular meshes.

We state the following conjecture:

Proposition 1 As h tends toward 0, the maximum pointwise error between the estimated laplacian described in Eq.5 and the real Laplace-Beltrami operator on the manifold tends toward 0 with a convergence speed in O(h

3 ), if one chooses parameter t in Θ(h 2 3 ).

Figure 3: The length of a linel l is drawn in orange as the result of the dot product between the estimated normal n and the elementary normal ne. It is an approximation of the projected arclength, drawn in blue.

Convergence results

We have implemented this digital Laplace-Beltrami operator with the DGtal library [DGt]. Input continuous smooth shapes are given as parameterized 2D closed curves , which are digitized using Gauss process (see an example in Fig. 4). Normal vectors along the digitized boundary are then estimated. These vectors are used to compute both the length between any two points and the measure of each point. Finally, we compare the real laplacian defined on the input curve with our estimated laplacian (from Eq.5) on the digitized curve.

We wish to confirm Prop.1. The choice of t is given by the following conditions: it must be as close as possible to 0, but not too close, otherwise we would fall back into a local laplacian, and convergence would be impossible. Therefore we choose t to be some h 2 3 for h the grid step. This choice of t leads to a convergence in O(h 2 3 ) as shown in Fig. 5.

Conclusion and discussion

We derived from [BSW08] a convolutional laplacian based on the heat equation. We showed its digitization and implementation on 2D digitized curves using Gauss process. We have checked Prop.1 by experimentation. We did not provide theoretical analysis, as it is still an ongoing work: the idea is to reformulate the theorem stated in [START_REF] Thibert | Properties of gauss digitized sets and digital surface integration[END_REF] and combine it with gaussian integration properties. Further work will include the implementation of this laplacian on 2D surfaces embedded in 3D as well as theoretical convergence proofs.

Références

Figure 1 :

 1 Figure 1: Implementation of geodesic distance computation described in [CWW13] using a discrete exterior calculus on digital surfaces. Isocontour in white on the left staired face of the cube shows metrics failure.

Figure 4 :

 4 Figure 4: An example of a 2D curve in blue and its digitization in black with h = 0.1. Pointels are black dots, linels are black line segments.

  into a local one c JFIG 2016.

Figure 5 :

 5 Figure 5: Convergence results using t = 23 on accelerated flower shape. Graphs in logarithm scale are to be read from right to left: as h tends toward 0, the error on laplacian estimation decreases with speed h 2 3 .

  [AFW06a] ARNOLD D. N., FALK R. S., WINTHER R.:Differential complexes and stability of finite element methods i. the de rham complex. In Compatible spatial discretizations. Springer, 2006, pp. 23-46. [AFW06b] ARNOLD D. N., FALK R. S., WINTHER R.: Finite element exterior calculus, homological techniques, and applications. Acta numerica. Vol. 15 (2006), 1-155. [BSW08] BELKIN M., SUN J., WANG Y.: Discrete laplace operator on meshed surfaces. In Proceedings of the 24th ACM Symposium on Computational Geometry, College Park, MD, USA, June 9-11, 2008 (2008), Teillaud M., (Ed.), ACM, pp. 278-287. [CLL14] COEURJOLLY D., LACHAUD J., LEVALLOIS J.: Multigrid convergent principal curvature estimators in digital geometry. Computer Vision and Image Understanding. Vol. 129 (2014), 27-41.

c journées JFIG 2016

c JFIG 2016.


Fabien Sixdenier
email: fabien.sixdenier@univ-lyon1.fr

Riccardo Scorretti
email: riccardo.scorretti@univ-lyon1.fr

Numerical model of static hysteresis taking into account temperature

Some soft magnetic materials are strongly dependent of the temperature, because of their low Curie temperature. In order to predict their behaviour in electrical devices, hysteresis models able to take into account the temperature are needed. The Vector Play static hysteresis model is an interesting hysteresis model that can be modified in order to take into account the thermal behaviour. In this paper, the temperature is taken into account in this last model through its parameters. The variations of some parameters with temperature are mainly issued from numerical interpolation and specific assumptions. Simulation results are compared to measurements and discussed.

Introduction

The well known Jiles-Atherthon [START_REF] Jiles | Ferromagnetic hysteresis[END_REF] and Preisach [START_REF] Preisach | Über die magnetische nachwirkung[END_REF] hysteresis models were widely used during the last decades to predict hysteresis loops or magnetic losses in many cases. Despite the huge number of users of these models, very few authors tried to take into account the temperature. For example [START_REF] Raghunathan | Theoretical model of dependence of hysteresis based on mean field theory[END_REF] used analytical empirical laws of variation of the Jiles-Atherthon model parameters versus temperature to predict hysteresis loops. In [START_REF] Messal | Temperature dependent extension of the jiles-atherton model: Study of the variation of microstructural hysteresis parameters[END_REF], authors compared several approaches (analytical [START_REF] Raghunathan | Theoretical model of dependence of hysteresis based on mean field theory[END_REF], numerical, and mixed) to identify each parameter variations. Concerning Preisach models, the authors of [START_REF] Stancu | Temperature and time-dependent preisach model for a stoner-wohlfarth particle system[END_REF] used a double well potential distribution. Both static and dynamic interactions for an ultra fine particulate system were taken into account in order to simulate various magnetization processes, such as: field-cooled, zero field-cooled, thermo-remanent magnetization, magnetic relaxation in applied field, and the major hysteresis loop as a function of temperature. The authors of [START_REF] Mitchler | Henkel plots in a temperature and time dependent preisach model[END_REF] assumed that the Preisach distribution function was a factorized product of a Gaussian coercive field distribution, and a Gaussian interaction field distribution both with some dispersion coefficients. Those last two approaches are rather complicated to implement in engineering softwares such as finite element or circuit-type (SPICE) ones.

The vector play models (also called energy-based vector hysteresis model) are phenomenological models of hysteresis first introduced by Bergqvist [START_REF] Bergqvist | Magnetic vector hysteresis model with dry friction-like pinning[END_REF] which have many desirable properties. In particular, it is intrinsically vectorial, and provides an accurate energy balance [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF]. It has been shown in [START_REF] François-Lavet | Vectorial incremental nonconservative consistent hysteresis model[END_REF] that the model accuracy was strongly improved by increasing the number of cells. However a relatively small number of internal variables was actually necessary to obtain an accurate representation of the complex response of the material. Therefore, a trade-off between accuracy and number of cells is always possible with this model. This is not the case for the Jiles Atherthon model for example.

The VPM has been tested and implemented in Finite Elements by many authors [START_REF] François-Lavet | Vectorial incremental nonconservative consistent hysteresis model[END_REF][START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF][START_REF] François-Lavet | An energy-based variational model of ferromagnetic hysteresis for finite element computations[END_REF] and recently been improved in the case of rotating fields [START_REF] Bergqvist | Temperature-dependent vector hysteresis model for permanent magnets[END_REF]. So as to be useful in practical engineering applications, models must consider the effect of temperature. The same authors also proposed a method [START_REF] Lin | Improved vector play model and parameter identification for magnetic hysteresis materials[END_REF] where they considered linear variations of saturation magnetization and coercive field to take into account the temperature for permanent magnets. In a previous paper [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF], variations of some parameters of the VPM were identified thanks to different papers in the literature [START_REF] Raghunathan | Theoretical model of dependence of hysteresis based on mean field theory[END_REF][START_REF] Weiss | La constante du champ moléculaire. equation d'état magnétique et calorimétrie[END_REF][START_REF] Gersdorf | Ferromagnetic properties of Fe and Ni in relation to their band structure[END_REF]. This last approach gave mixed results. First, the purely magnetic behaviour gave not so good results for very low applied fields (not enough accurate description in this particular region). Second, the VPM showed no difficulties to adapt itself to the temperature dependence, by allowing some parameters to change their values. The thermal behaviour was qualitatively well retrieved, loops expanded or shrank themselves in function of temperature showing the ability of the build model to qualitatively retrieve the Hopkinson effect.

In order to improve the accuracy of the simulations of the VPM, we present in this work recent developments on the dependence of VPM on temperature. In our approach, no analytical (empirical) laws are used in order to predict the thermal behaviour. Secondarily, a novel approach for identifying VPM's parameters is proposed.

The paper is organized as follows: the implemented VPM and the main assumptions made to take into account the temperature are presented (section 2). The identification of VPM's parameters is analysed in section 3. In section 4, the results of the identification algorithm are presented on a virtual and a real material. Discussion of obtained results and perspectives conclude the paper (section 5).

The vector play model 2.1 Purely magnetic behaviour

Hysteresis in magnetic materials can be seen as defects at microscopical scale, which pin Bloch walls and prevent them from moving under the effect of an external magnetic field. When the magnetic energy becomes higher than the energy barrier due to this pinning force, Bloch walls suddenly rearrange and the energy corresponding to the jump is dissipated [START_REF] Bergqvist | Magnetic vector hysteresis model with dry friction-like pinning[END_REF].

Basing on an analogy between pinning force and mechanical friction force, Henrotte et al. proposed a phenomenological (macroscopic) model which is consistent with laws of thermodynamics [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF]. In this model, the applied magnetic field h is decomposed in a reversible h re and irreversible h ir part:

h = h re + h ir
The magnetic field h is the input of the model. The reversible magnetic field h re is computed as the weighted sum of contributions h k re of a number N of "cells" (in practice h ir is never computed explicitly). Each cell is analogous to a mechanical system composed of a spring and a friction slider:

h k re = N k=1 ω k h k re (1)
In this sum, ω k ≥ 0 are the weights associated to cells, and must verify the condition:

N k=1 ω k = 1.
The terms h k re are the internal state of each cell. They are upgraded as a function of the applied magnetic field h and of their current value h k re0 according to:

h k re = h k re0 if ||h -h k re0 || < χ k h -χ k • h-h k re0 ||h-h k re0 || otherwise (2)
where χ k is a parameter linked to the coercivity [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF]. The relationship between h, h re and h i is depicted in figure 1. The magnetization M is computed as a function of the reversible magnetic field as:

M = M an (||h re ||) • h re ||h re || (3)
where M an ( |h re ||) is a scalar anhysteretic magnetization function. Finally, the magnetic flux density is computed as:

b = µ 0 (M + h) (4)
One of the advantages of this model is that it is intrinsically vectorial. Moreover, the number of cells is not fixed a priori, and can eventually be adjusted in such a way to fulfill accuracy requirements.

Temperature dependent extension

All equations and principles of the model seen before remain valid, only some parameters will have their value changed. As seen in the previous section, the VPM is mainly governed by the anhysteretic function M an and by the χ k values. For real materials, analytical anhysteretic functions (see [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF] for possible analytical functions) always present more or less error on one or several regions (saturation, saturation knee, low field region) compared to the measured anhysteretic curve. Moreover, analytical functions also need an identification algorithm of physical (best case) or empirical (worst case) parameters. As the temperature will change, all anhysteretic function parameters will change as well. If some physical parameters variations with temperature are well described in the literature (example : saturation magnetization M s (T )), some empirical parameters are not. For these reasons, we chose to use a specific interpolation table in order to evaluate M an at a given temperature.

For anhysteretic magnetization M an , several major loops are measured at different temperatures that were performed on a ring sample wound with high temperature enameled wires (placed in an thermal enclosure) following the fluxmetric method with respect to the IEC 60404-4 standard [START_REF] Fausto | Characterization and Measurement of Magnetic Materials[END_REF]. More precisely, the major loop is measured at several temperatures, and the anhysteretic magnetization M an (h, T ) is computed by considering the middle point of horizontal lines which connect ascending and descending branches of the (major) loop. The anhysteretic magnetization function M an (h, T ) for a generic temperature T is obtained by interpolation. The figure 2 (left) shows the built 2D anhysteretic magnetization M an versus the temperature T and of the applied field h. These data are measured on the same material than the one presented in section 4. As for parameters ω k , the same property (

N k=1 ω k = 1
) must hold at any temperature, we assume that in fact they don't depend on the temperature. However parameters χ k are related to a local pinning force or a local coercive field; therefore it is assumed here that the variations of all χ k with temperature are the same as the macroscopic coercive field H c , which is already known with the measurements of major loops at different temperatures (figure 2 (right)). Now, let us suppose that all χ k have been identified at room temperature (RT). All χ k at any temperature T within the range can be rescaled by applying [START_REF] Stancu | Temperature and time-dependent preisach model for a stoner-wohlfarth particle system[END_REF].

χ k (T ) = H c (T ) H c (RT ) χ k (RT ) (5)
In order to verify the assumptions made here, the (ω k , χ k ) parameters have to be identified so as to that the VPM give accurate results for any kind of signals (low/high levels, single/multiple frequency(ies)).

Parameters identification

The VPM is fully characterized by the anhysteretic magnetization function M an (h), and the N couples of parameters (ω k , χ k ) which represent each of the N cells which compose the model. The number N of cells is a priori unknown.

The parameters ω k and χ k must verify the following constraints:

ω k ≥ 0 ; N k=1 ω k = 1 ( 6
)
χ 1 = 0 ; χ k ≤ χ k+1 (7)
In [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF], it is proposed to determine these parameters basing on the experimental curve H c (H p ), where H c is the coercitive field corresponding to a minor symmetric closed loop, for which the peak magnetic field is H p . The same approach is followed in [START_REF] François-Lavet | An energy-based variational model of ferromagnetic hysteresis for finite element computations[END_REF] and [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF]. However, this approach requires a specific set of measurements (i.e. symmetric closed loops at many peak levels), which have to be performed purposely. Conversely, it would be interesting to explore the possibility of being able to identify model parameters from a generic set of measurements. In this work we followed this second approach. We stress that the term "generic" does not mean that any dataset will suit, but rather that the excitation field h(t) does not need to be a specific signal. This last point will be discussed later.

Anhysteretic magnetization

The anhysteretic magnetization function M an (h, T ) is obtained from experimental measurements taken at different temperatures (see previous section).

Naive formulations for the identification problem

Therefore the only parameters which have to be identified are ω k and χ k (at a given temperature T ) with k spanning from 1 and the (unknown) number of cells N . According to [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF], these parameters have been fitted for a single temperature T 0 = 298 K; their variation with temperature is as described in section 2. For sake of simplicity, let's define the vectors Y = (ω 1 . . . ω N ) and X = (χ 1 . . . χ N ). For a given number of cells N , and a given dataset of experimental measurement b(h), one possibility would be to solve the constrained optimization problem:

     (ω k , χ k ) = arg min Y,X ||b(h) -b s (h)|| 2 2 ω k ≥ 0 ; N k=1 ω k = 1 χ 1 = 0 ; χ k ≤ χ k+1 (8)
where b s (h) is the magnetic flux density computed by the model, which clearly depends on the set of parameters (ω k , χ k : k = 1 . . . N ). In fact, the quadratic term ||b(h) -b s (h)|| 2 could be replaced by any meaningfull cost function, for instance:

     (ω k , χ k ) = arg min Y,X ||H c (H p ) -H c,s (H p )|| 2 2 ω k ≥ 0 ; N k=1 ω k = 1 χ 1 = 0 ; χ k ≤ χ k+1 (9)
where H c,s (H p ) is the simulated coercive field, which can be easily evaluated [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF] as:

H c,s (H p ) = m(Hp) k=1 ω k χ k m(Hp) k=1
ω k [START_REF] François-Lavet | An energy-based variational model of ferromagnetic hysteresis for finite element computations[END_REF] where the index m(H p ) ∈ [1; N ] is the largest integer such that χ k ≤ H p ∀k ≤ m(H p ).

Improved formulation

All these "naive" formulations are prone to several problems. One of them is that the optimization algorithm can stagnate in a local minimum. Before introducing better formulations, one can observe that a single cell (ω, χ) can be replaced by two cells (ω 1 , χ 1 ) and (ω 2 , χ 2 ), provided that: ω = ω 1 + ω 2 and χ = χ 1 = χ 2 . With more generality, a single cell (ω, χ) can be replaced by a number p of cells (ω i , χ) in the following way:

(ω, χ) ⇔ (ω 1 , χ) . . . (ω p , χ) provided that:

p i=1 ω i = ω (11)
Keeping this observation in mind, we can reformulate the problem (8) by considering a number N of cells (ω, χ k ), all of which share the same value of ω = 1 N : in this way, the constraints (6) are fulfilled by construction:

     (ω k , χ k ) = arg min X ||b(h) -b s (h)|| 2 2 ω k ≡ 1/N χ 1 = 0 ; χ k ≤ χ k+1 (12)
The formulation [START_REF] Lin | Improved vector play model and parameter identification for magnetic hysteresis materials[END_REF] can be further simplified by taking as unknown the difference ∆χ k = χ k+1 -χ k instead of χ k , which will be constrained to be nonnegative. In this way we obtain the following (nonlinear) non-negative least squares formulation:

∆χ k = arg min ∆X ||b(h) -b s (h)|| 2 2 ∆χ k ≥ 0 ( 13
)
where ∆X = (∆χ 1 . . . ∆χN -1). Observe that formulation (8) has a number 2N -1 unknowns, whereas this formulation has only N -1 unknowns. However, due to the fact that in this formulation ω k takes discrete values, in principle a finer discretization (i.e. a higher value of N ) would be required so as to achieve comparable accuracy.

The most important improvement of formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] with respect of (8) or ( 9) is that constraints are much more relaxed. Moreover, due to the equivalence property [START_REF] Bergqvist | Temperature-dependent vector hysteresis model for permanent magnets[END_REF], the fact that in (8) ω k are free parameters may introduce many local minima in the cost function. This is not the case with [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] because ω k ≡ 1/N is imposed. Therefore the optimization problem should be easier to solve numerically.

The price to pay for this is that a higher number N of cells has to be considered. This point is a huge drawback in view of Finite Element implementation. Therefore we devised an algorithm for reducing the number of cell. Namely, our strategy is to solve [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] with a high number of cells, and then simplify the obtained VPM by packing some cells together.

Packing cells algorithm

First of all, observe that the number of independent cells provided by formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] corresponds to the number of non-zero values of ∆X. In fact, two (or more) consecutive cells corresponding to ∆χ k = 0 can be packed together into a single cell by using the same argument [START_REF] Bergqvist | Temperature-dependent vector hysteresis model for permanent magnets[END_REF] developed beforehand, that is:

∆χ k = . . . = ∆χ k+p-1 = 0 ⇔          (1/N, χ k ) (1/N, χ k+1 ) . . . (1/N, χ k+p ) =          (1/N, χ k ) (1/N, χ k ) . . . (1/N, χ k ) -→ ( p N , χ k ) (14)
The main idea of the algorithm is to consider the case of a VPM composed by an unknown, small number N T of independent cells, which is identified by solving the problem [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF]. In the ideal case (no measurement noise, exact arithmetic, etc.) the solution (∆χ k ) should be composed of exactly N T non-zero terms, whatever the discretization N . However in practice this is rarely the case. So, we are expecting that the obtained solution will be composed of at best N T non-zero terms, and other small but non-zero terms.

In order to recover the true solution, the idea is to search the smallest nonzero term, and to aggregate it with its bigger neighbour. Namely, assume that ∆χ k is the smallest non-zero term of ∆X, and that ∆χ k+1 > ∆χ k-1 ; then the algorithm will execute:

∆χ k+1 ← ∆χ k+1 + ∆χ k ∆χ k ← 0
This simple iteration is repeated until it is possible to simplify the VPM, that is until two non-zero, not too small, adjacent terms can be found in the vector ∆X. At each iteration a term of ∆X is set to 0, so that this algorithm will necessarily terminate. The pseudocode of the algorithm is listed hereafter:

1: procedure simplifyModel(∆X[1 . . . N -1]) 2: repeat 3: find I = {i ∈ [1, N -1] : ∆χ i = 0 ∧ max(∆χ i+1 , ∆χ i-1 ) > 0} 4: if I = ∅ then 5: i ← arg min i∈I (∆χ i ) 6: j ← arg max j∈{i-1,i+1}
(∆χ j ) 7:

∆χ j ← ∆χ j + ∆χ i 8:

∆χ i ← 0 9:
end if Table 1: Values of χ k (A/m) obtained from simulations (N = 10) k eq. ( 13) eq. ( 13), simplified Table 2: Values of χ k (A/m) obtained for a M nZn N30 ferrite (N = 10) k eq. ( 13) eq. ( 13), simplified simplification algorithm is very effective in reducing the number of cells of the VPM, yet ensuring a good agreement with data. Computations performed with higher values of N provided quite similar results (figure 5 left). We compared our approach with the one proposed in [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF]. We found that the identified VPMs are compatible with the curve H c (H p ) for low values of H p (namely H p < 60 A/m), whereas discrepancy is observed for higher values (figure 5

right).

We tried to identify the VPM by using a major cycle only, but we obtained a very poor result (not shown). We believe that this is due to the fact that cells with lowest values of χ k , which come to play with small loops, are the hardest to identify. When using large loops for identification, these cells come to play only when the field is reversed -that is for high values of the flux density, where the sensibility of measurements is low. Therefore, it can be argued that it is necessary to include small loops in the dataset used for identification.

Various tests of the VPM with h(t) and T

We tested the VPM with several kind of signals h and temperatures T . The set of parameters at room temperature are those mentioned in table 2, with formulation (13) simplified and the anhysteretic magnetization functions M an (h, T ) and χ(T ) that are used are those plotted in Fig 2 . Figure 6 (left) shows measured and simulated major and minor symmetric loops at room temperature showing that the VPM is able to reproduce symmetric loops at room temperature with different levels of applied field. Figure 6 (right) shows that the VPM can simulate major loops at different temperatures thanks to the method mentioned in section 2.

As the VPM is normally able to predict loops with any kind of h(t) signals, we tested it with the applied field defined by [START_REF] Gersdorf | Ferromagnetic properties of Fe and Ni in relation to their band structure[END_REF]:

h(t) = H max (sin (ωt) + sin (3ωt)) (15)
Simulated and measured loops are compared on figure 7 for two levels of H max . Simulated and measured results are in very good agreement showing that the set of identified parameters combined to the method defined in section 2 is Last but not least, we wanted to know if the VPM combined to its dependent extension was able to retrieve the Hopkinson effect [START_REF] Hopkinson | Magnetic and other physical properties of iron at a high temperature[END_REF]. Loops are simulated at different temperatures for two applied field levels (h imposed). The tips of each loop is used to determine the per unit magnetic permeability µ r . Figure 8 (left) show µ r variations with the temperature T for two levels of applied field (H max = 8 A/m and H max = 40 A/m). When a low field (H max = 8 A/m) is applied to the material, a peak of magnetic permeability µ r is observed when T is increasing just before the Curie temperature, then it drops to µ r ≈ 1, at T c (Curie temperature). Figure 8

H max =400 A/m (m) H max =400 A/m (s) H max =55 A/m (m) H max =55 A/m (s) H max =30 A/m (m) H max =30 A/m (s) H max =15 A/m (m) H max =15 A/m (s) H max =8 A/m (m) H max =8 A/m (
H (A/m) 223 K (m) 248 K (s) 298 K (m) 298 K (s) 348 K (m) 348 K (s) 398 K (m) 398 K (s) 423 K (m) 423 K (s) -0.1 0 0.1 -10 -5 0 5 10 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4

Discussion

The main purpose of this work is to describe how Vector Play Model (VPM) can effectively take into account the effect of temperature. In particular, it has been shown that VPM identified at room temperature can be used to model the behaviour of magnetic materials at any temperature, provided that few additional measurements have been performed in the range of temperature.

To this aim, we developed formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] to identify the parameters of the VPM. Conversely to existing works [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF][START_REF] François-Lavet | An energy-based variational model of ferromagnetic hysteresis for finite element computations[END_REF][START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF], this formulation has as unknowns only the parameters (χ k ), whereas ω k ≡ 1/N is fixed a priori.

As for the measurements required to identify VPMs, formulation (13) requires no particular signal: the only requirement is that datasets have to be "rich enough" for activating all cells in such a way that they can be discriminated. For instance, symmetric minor cycles will suit for this purpose, whereas the major cycle alone will not -the reason being that cells corresponding to small values of χ k are activated for high values of the field, where the sensitivity is the poorest: so they can be hardly identified.

We compared our approach with the approach proposed in [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF]. We found that, according to our measurements, formulation (13) provides more accurate result, which is compatible with the H c (H p ) curve for low values of H p . However, it has to be said that in our work VPMs have been identified by using a cycle, the peak excitation field of which is 45 A/m: therefore it could be reasonable that some discrepancy may appear for higher values of H p . This point merits to be more deeply investigated.

It has to be remarked that this formulation is very stable with respect of the number of cells N . The approach followed in this work is therefore to identify VPM with a high number of cells, and then simplify it by using the proposed (or others) algorithm, which has demonstrated to be very effective in reducing the complexity of the model. Moreover, this algorithm can be further improved: for instance, in the exemple discussed above it is clear that the last two cells could be packed together to obtain a 3-cells VPM: By the way, this approach does not require the a priori knowledge of the number of cells (which is needed with formulation (8)), or -said differently -it is possible to devise algorithms which determine the minimum number of cells which is required to achieve an imposed accuracy.

It is shown (figures 6 and 7) that the experimentally measured loops were very well retrieved by the VPM. This fact suggests that (at least for this material) the assumptions made at section 2 regarding the thermal behaviour are correct. It also points out that the χ k identification protocol is a key element in order to have a good accuracy in a large range.

Concerning the Hopkinson effect, the VPM behaves qualitatively well due, we think to the following reasons. For low field levels, at low temperature only the smallestχ k are activated (moving). When T increase, all χ k decrease and new cells can then be activated so the relative permeability µ r increase until the material saturates. The peak of µ r can then be seen as new cells that are activated in the material as T increase. When the material is saturated, µ r decrease because the saturation magnetization decrease very rapidly as T is close to the Curie temperature (this is due to the dependence of M an (h re , T ) on T ).

For higher field levels, all "small" χ k cells are activated even at low temperature. The small increase of µ r is, in this case, due to the small reduction of the large χ k cells which tend to make the cells a bit more sensitive to the field. But no new cells are activated as T increase. The decrease of µ r is also due to the saturation, but arrives at lower temperature because the field level is higher.

Conclusion and perspectives

In this work we propose modifications to existing VPM so as to take into account the effect of temperature, as well as a new formulation for identifying VPM from experiments. Conversely to existing approaches [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF][START_REF] Krah | Numerical optimization of a hysteresis model[END_REF], our approach does not require a specific kind of measurement. Proposed formulation has been shown to be stable, and by using a simplification algorithm can provide accurate VPMs with a low number of cells. Some discrepancy appeared between our results and the approach proposed by [START_REF] Henrotte | An energy-based vector hysteresis model for ferromagnetic materials[END_REF], which require some more work to be elucidated. The VPM associated with a 2D Lookuptable for M an and a rescaling of the χ k thanks to the macroscopic coercive field H c is now able to predict any loop for any kind of signals h(t) for any temperature within the range. VPM has also shown its ability to retrieve qualitatively the Hopkinson effect.

Up to now, the proposed approach has been tested successfully with a single material: future works will be devoted to more extensive experimentation with different materials, so as to assess (or not) the effectiveness of our model.

Possible other extensions of the VPM model, following the same approach, may be envisaged. For example, simulate the influence of a mechanical stress on the magnetic properties by modifying M an and χ k following the same method and assumptions than here (lookup tables and rescaling).

Another possible extension would be to simulate dynamic hysteresis. Two separate ways can be imagined. The first one, would be to decompose the dynamic applied field in several fields such as Bertotti's theory [START_REF] Bertotti | General properties of power losses in soft ferromagnetic materials[END_REF] and to use the VPM to compute the static part of the field H s . To do that, the VPM has to be inversed i.e calculate H s (B). Another way, would be to attribute a mass to the different cells in order to give inertia to the cells.

1 Figure 1 :

 11 Figure 1: Vector diagram of h, h re and h i for a single cell (left), and for N = 3 cells (right).

2 .Figure 2 :

 22 Figure 2: Left: Anhysteretic magnetization M an versus the temperature T and of the applied field H . Right: Coercive field H c versus temperature T . Measurements were made on a M nZn N30 ferrite material (see section 4).

Figure 3 :

 3 Figure 3: Left: simulated noisy measurement and fitted VPM by using formulation 13 with N = 10. Right: simulation performed on a validation dataset with the ordinary (red) and simplified (yellow) model; in the inset is shown a zoom of a minor asymmetrical loop.

Figure 4 :

 4 Figure 4: Left: measurement used to fit VPMs with N = 10. Right: validation with a symetric minor loop. Red and yellow curves are the ouput generated respectively by the ordinary (red) and simplified (yellow) VPM.

Figure 5 :

 5 Figure 5: Left: estimated values of (χ k ) for VPMs with N = 10, 20 and 40 cells. Right: curve H c (H p ) computed from VPMs by using (10) (lines) and experimental measurements (+).

Figure 6 :

 6 Figure 6: Left : major and minor symmetric loops at room temperature. Right : Major loops for different temperatures. Full lines : measured values ; dotted lines : simulated values

Figure 7 :

 7 Figure 7: Left : loops with third harmonic (H max = 40 A/m) at different temperatures. Right : loops with third harmonic (H max = 80 A/m) at different temperatures. Full lines : measured values ; dotted lines : simulated values

Figure 8 :

 8 Figure 8: Left : relative permeability µ r vs T . Right : minor loops at different temperatures showing the Hopkinson effect.

  (0.1, 0) -(0.3, 5.064) -(0.5, 14.938) -(0.1, 15.454) (0.6, 15.024)

10:

until I = ∅ 11: end procedure 4 Results

Fake material

Before using formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] to identify a real material, we performed preliminary tests by using simulated measurements and validation datasets. The fake material we simulated has an anhysteretic magnetization described by a Langevin function: The data generated by using this model have been fitted with formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF]. So as to avoid inverse crime [START_REF] Wirgin | The inverse crime[END_REF], in the generation of this synthetic data, we simulated a Gaussian white noise η ∼ N (0, 5 2 mT 2 ), and we introduced a bias of 1% in parameters M s and h 0 . The minimization algorithm used is the interior point algorithm [START_REF] Byrd | An interior point algorithm for large-scale nonlinear programming[END_REF], as implemented in the function fmincon of Matlab R . All simulations have been repeated 100 times so as to evaluate confidence intervals by a Monte Carlo method [START_REF] Aster | Parameter estimation and inverse problems[END_REF]. The result of one of these is depicted in figure 3. The values of (χ k ) obtained for N = 10 are listed in One observes that all formulations provide results which are quite close to the reference VPM (figure 3). VPMs with different numbers of cells can be compared by plotting χ k as a function of the variable Ω k = k i=1 ω i . We performed simulations with N spanning from 10 to 50: we found that the obtained results are quite stable (not shown).

The simplification algorithm is very effective in reducing the number of cells. For engineering purposes, one could mention that the computation of losses also provide accurate results up to 5% for all fitted models.

Real material

We identified different VPMs corresponding to an Manganese Zinc N30 ferrite material, which has been selected for its low value of T c < 433 K. The experimental measurements used for identification and validation are depicted in figure 4. For identification we used the formulation (13), with N spanning from 10 to 50. Datasets used for identification and validation have been measured at 298 K. The obtained values for (χ k ) are reported in table 2 for N = 10.

One observes that the formulation [START_REF] Sixdenier | Temperature-dependent extension of a static hysteresis model[END_REF] provided results in good agreement with experimental data (comparison of loops). Like in previous simulations, the


Abdallah Saffidine
email: abdallah.saffidine@dauphine.fr

Tristan Cazenave
email: cazenave@lamsade.dauphine.fr

Developments on Product Propagation

Introduction

Product Propagation (PP) is a way to backup probabilistic information in a two-player game tree search [START_REF] Slagle | Experiments with a multipurpose, theoremproving heuristic program[END_REF]. It has been advocated as an alternative to minimaxing that does not exhibit the minimax pathology [START_REF] Pearl | On the nature of pathology in game searching[END_REF][START_REF] Baum | Best play for imperfect players and game tree search[END_REF][START_REF] Baum | A Bayesian approach to relevance in game playing[END_REF][START_REF] Horacek | Towards understanding conceptual differences between minimaxing and product-propagation[END_REF][START_REF] Horacek | An analysis of decision quality of minimaxing vs. product propagation[END_REF].

PP was recently proposed as an algorithm to solve games, combining ideas from Proof Number Search (PNS) and probabilistic reasoning [START_REF] Stern | Learning to solve game trees[END_REF]. In Stern's paper, PP was found to be about as performant as PNS for capturing go problems.

We conduct a more extensive study of PP, comparing it to various other paradigmatic solving algorithms and improving its memory consumption and its solving time. Doing so, we hope to establish that PP is an important algorithm for solving games that the game search practician should know about. Indeed, we exhibit multiple domains in which PP performs better than the other tested game solving algorithms.

The baseline game tree search algorithms that we use to establish PP's value are Monte Carlo Tree Search (MCTS) Solver [START_REF] Mark | Monte-Carlo tree search solver[END_REF] which was recently used to solve the game of havannah on size 4 [START_REF] Ewalds | Playing and solving Havannah[END_REF]; PNS [START_REF] Victor | Proof-Number Search[END_REF][START_REF] Van Den Herik | Proof-Number Search and its variants[END_REF][START_REF] Kishimoto | Game-tree search using proof numbers: The first twenty years[END_REF]; and αβ [START_REF] Russell | Artificial Intelligence -A Modern Approach[END_REF].

The next section deals with algorithms solving two-player games. The third section is about the Product Propagation algorithm. The fourth section details experimental results and the last section concludes.

Solving two-player games

We assume a deterministic zero-sum two-player game with two outcomes and sequential moves. Game tree search algorithms have been proposed to address games with multiple outcomes [START_REF] Cazenave | Score bounded Monte-Carlo tree search[END_REF][START_REF] Saffidine | Multiple-outcome proof number search[END_REF], multi-player games [START_REF] Sturtevant | A comparison of algorithms for multi-player games[END_REF][START_REF] Maarten | Best reply search for multiplayer games[END_REF][START_REF] Nijssen | An overview of search techniques in multi-player games[END_REF], non-deterministic games [START_REF] Hauk | Rediscovering *-minimax search[END_REF], and games with simultaneous moves [START_REF] Saffidine | Alpha-beta pruning for games with simultaneous moves[END_REF].

A generic best-first search framework is presented in Algorithm 1. To instantiate this framework, one needs to specify a type of information to associate to nodes of the explored tree as well as functions to manipulate this type: info-term, init-leaf, select-child, update.

PNS is a best-first search algorithm which expands the explored game tree in the direction of most proving nodes, that is, parts of the tree which seem easier to prove or disprove.

MCTS Solver is also a best-first search algorithm which can be cast in the mentioned framework. In MCTS [START_REF] Browne | A survey of Monte Carlo tree search methods[END_REF], the information associated to nodes is in the form of sampling statistics, and bandit-based formulas [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF] are used to guide the search. The sampling performed at a leaf node in MCTS can take the form of games played randomly until a terminal position, but it can also be the value of a heuristical evaluation function after a few random moves [START_REF] Lorentz | Amazons discover Monte-Carlo[END_REF][START_REF] Mark | Monte Carlo tree search in lines of action[END_REF]. We denote the latter variant as MCTS-E.

Probability Search

The Product Propagation Algorithm

Product Propagation (PP) is a recently proposed algorithm to solve perfect information two-player two-outcome games based on an analogy with probabilities [START_REF] Stern | Learning to solve game trees[END_REF].

In PP, each node n is associated to a single number PPN(n) (the probability propagation number for n) such that PPN(n) ∈ [0, 1]. The PPN of a leaf corresponding to a Max win is 1 and the PPN of a Max loss is 0. PPN(n) can intuitively be understood as the likelihood of n being a Max win given the partially explored game tree. With this interpretation in mind, natural update rules can be proposed. If n is an internal Min node, then it is a win for Max if and only if all children are win for Max themselves. Thus, the probability that n is win is the joint probability that all children are win. If we assume all children are independent, then we obtain that the PPN of n is the product of the PPN of the children for Min nodes. A similar line of reasoning leads to the formula for Max nodes. To define the PPN of a non-terminal leaf l, the simplest is to assume no information is available and initiate PPN(l) to 1 2 . These principles allow to induce a PPN for every explored node in the game tree and are summed up in Table 1.

Note that this explanation is just a loose interpretation of PPN(n) and not a formal justification. Indeed, the independence assumption does not hold in practice, and in concrete games n is either a win or a loss for Max but it is not a random event. Still, the independence assumption is used because it is simple and the algorithm works well even though the assumption is usually wrong.

To be able to use the generic best first search framework, we still need to specify which leaf of the tree is to be expanded. The most straightforward ap-bfs(state q, player m) r ← new node with label m r.info

← init-leaf(r) n ← r while r is not solved do while n is not a leaf do n ← select-child(n) extend(n) n ← backpropagate(n) return r extend(node n) switch on the label of n do case terminal n.info ← info-term(n) case max foreach q in {q , q a -→ q } do n ← new node with label min n .info ← init-leaf(n ) Add n as a child of n case min foreach q in {q , q a -→ q } do n ← new node with label max n .info ← init-leaf(n ) Add n as a child of n backpropagate(node n) new_info ← update(n) if new_info = n.info ∨ n = r then return n else n.info ← new_info return backpropagate(n.parent)
Algorithm 1: Pseudo-code for a best-first search algorithm.

Table 1. Initial values for leaf and internal nodes in PP. C denote the set of children.

Node label PPN info-term

Max wins

1 Max loses 0 init-leaf 1 2 update Max 1 -C (1 -PPN) Min C PPN
proach is to select the child maximizing PPN when at a Max node, and to select the child minimizing PPN when at a Min node, as shown in Table 2.

Practical improvements

The mobility heuristic provides a better initialization for non-terminal leaves. Instead of setting PPN to1 /2 as described in Table 1, we use an initial value that depends on the number of legal moves and on the type of node. Let c be the number of legal moves at a leaf, the PPN of which we want to initialize. If the leaf is a Max -node, then we set

PPN = 1 -1 /2 c . If the leaf is a Min-node, then we set PPN = 1 /2 c .
In the description of best first search algorithms given in Algorithm 1, we see that new nodes are added to the memory after each iteration of the main loop in bfs. Thus, if the init-leaf procedure is very fast then the resulting algorithm will fill the memory very quickly. Earlier work on PNS provides inspiration to address this problem [START_REF] Kishimoto | Game-tree search using proof numbers: The first twenty years[END_REF]. For instance, Kishimoto proposed to turn PP into a depth-first search algorithm with a technique similar to the used in dfpn. 1 Alternatively, it is possible to adapt the PN 2 ideas to develop a PP 2 algorithm. In PP 2 , instead of initializing directly a non-terminal leaf, we call the PP algorithm on the position corresponding to that leaf with a bound on the number of nodes. The bound on the number of nodes allowed in the sub-search is set to the number of nodes that have been created so far in the main search. After a sub-search is over, the children of the root of that search are added to the tree of the main search. Thus, the PPN associated these newly added nodes is based on information gathered in the sub-search, rather than based only on an initialization heuristic.

Experimental Results

While the performance of PP as a solver has matched that of PNS in go [START_REF] Stern | Learning to solve game trees[END_REF], it has proven to be disappointing in shogi. 1 We now exhibit several domains where the PP search paradigm outperforms more classical algorithms.

In the following sets of experiments, we do not use any domain specific knowledge. We are aware that the use of such techniques would improve the solving ability of all our programs. Nevertheless, we believe that showing that a generic and non-optimized implementation of PP performs better than generic and nonoptimized implementations of PNS, MCTS, or αβ in a variety of domains provides good reason to think that the ideas underlying PP are of importance in game solving.

We have described a mobility heuristic for PP variants in Section 3.2. We also use the classical mobility heuristic for PNS variants. That is, if c is the number of legal moves at a non-terminal leaf to be initialized, then instead of setting the proof and disproof numbers to 1 and 1 respectively, we set them to 1 and c if the leaf is a max -node or to c and 1 if the leaf is a min-node.

All variants of PNS, PP, and MCTS were implemented with the best-first scheme described in Section 2. For PN 2 and PP 2 , only the number of nodes in the main search is displayed.

The game of y

The game of y was discovered independently by Claude Shannon in the 50s, and in 1970 by Schensted and Titus [26]. It is played on a triangular board with a hexagonal paving. Players take turns adding one stone of their color on empty cells of the board. A player wins when they succeed in connecting all three edges with a single connected group of stones of their color. Just as hex, y enjoys the no-draw property.

The current best evaluation function for y is the reduction evaluation function [START_REF] Van Rijswijck | Search and evaluation in Hex[END_REF]. This evaluation function naturally takes values in [0, 1] with 0 (resp. 1) corresponding to a Min (resp. Max ) win.

PNS with the mobility initialization could not solve any position in less than 3 minutes in a preliminary set of about 50 positions. As a result we did not include this solver in our experiment with a larger set of positions. The experiments on y was carried out as follows. We generated 77,012 opening positions on a board of size 6. We then ran PP using the reduction evaluation function, MCTS using playouts with a random policy, and a variant of MCTS using the same reduction evaluation instead of random playouts (MCTS-E). For each solver, we recorded the total number of positions solved within 60 seconds. Then, for each solving algorithm, we computed the number of positions among those 77,012 which were solved faster by this solver than by the two other solver, as well as the number of positions which needed fewer iterations of the algorithm to be solved. The results are presented in Table 3.

We see that the PP algorithms was able to solve the highest number of positions, 77,010 positions out of 77,012 could be solved within 60 seconds. We also note that for a very large proportion of positions (68,477), PP is the fastest algorithm. However, MCTS needs fewer iterations than the other two algorithms on 35,444 positions. A possible interpretation of these results is that although iterations of MCTS are a bit more informative than iterations of PP, they take much longer. As a result, PP is better suited to situations where time is the most important constraint, while MCTS is more appropriate when memory efficiency is a bottleneck. Note that if we discard MCTS-E results, then 72,830 positions are Figure 1 displays some of these results graphically. We sampled about 150 positions of various difficulty from the set of 77,012 y positions, and plotted the time needed to solve such positions by each algorithm against the time needed by PP. We see that positions that are easy for PP are likely to be easy for both MCTS solvers, while positions hard for PP are likely to be hard for both other solvers as well.

domineering

domineering is played on a rectangular board. The first player places a vertical 2 × 1 rectangle anywhere on the board. The second player places an horizontal 2 × 1 rectangle, and the games continues like that until a player has no legal moves. The first player that has no legal moves has lost.

domineering has already been studied in previous work by game search specialists as well as combinatorial game theorists [START_REF] Breuker | Solving 8×8 domineering[END_REF][START_REF] Michael Lachmann | Who wins domineering on rectangular boards[END_REF]. 2 While these papers focusing on domineering obtain solution for relatively large boards, we have kept ourselves to a naive implementation of both the game rules and the algorithms. In particular, we do not perform any symmetry detection nor make use of combinatorial game theory techniques such as decomposition into subgames.

We presents results for the following algorithms: αβ, PNS with Transpositions (PNT) [START_REF] Schijf | Proof-number search and transpositions[END_REF], PN 2 [START_REF] Michel | Memory versus Search in Games[END_REF], PP, PP with Transpositions (PPT) and PP 2 . The PNS algorithm could not find a single solution within 10 7 node expansion when transpositions where not detected and it is thus left out.

For PNS variants the standard mobility heuristic is used to compute the proof numbers and the disproof numbers at non solved leaves. For PP variants, we used the mobility heuristic as described in Section 3.2.

Tables 4 and5 give the number of nodes and times for different algorithms solving domineering. αβ is enhanced with transposition tables, killer moves, the history heuristic and an evaluation function. We can see that on the smallest 5 × 6 board that PPT gives the best results. On the larger 6 × 6 board PPT is the best algorithm by far. On the largest 7 × 6 board, several algorithms run out of memory, and the best algorithm remains PPT which outperforms both αβ and PN 2 . In their paper, Breuker et al, have shown that the use of transposition tables and symmetries increased significantly the performance of their αβ implementation [START_REF] Breuker | Solving 8×8 domineering[END_REF]. While, our proof-of-concept implementation does not take advantage of symmetries, our results show that transpositions are of great importance in the PP paradigm as well. For standard board sizes such as 4 × 4 or 5 × 4, αβ gives the best results among the algorithms we study in this paper. We have noticed that for N × 1 boards for N ≥ 18, PPT becomes competitive. Results for a few board sizes are given in Table 6 for the number of nodes and in Table 7 for the times.

Conclusion

In this paper, we have presented how to use Product Propagation (PP) in order to solve abstract two-player games. We extended PP so as to handle transpositions and to reduce memory consumption with the PP 2 algorithm. For two of the games that have been tested (i.e., y, domineering), we found that our extensions of PP are able to better solve games than the other solving algorithms. For nogo, PP variants outperform PNS variants on all tested sizes, and PP does better than αβ on some sizes but αβ is better on standard sizes. Being a best-first search algorithm, PP is quite related to PNS and MCTS, as such, it seems natural to try and adapt ideas that proved successful for these algorithms to the Product Propagation paradigm. For instance, while PNS and PP are originally designed for two-outcome games, future work could adapt the ideas underlying Multiple-Outcome PNS [START_REF] Saffidine | Multiple-outcome proof number search[END_REF] to turn PP into an algorithm addressing more general games. Adapting more elaborate schemes for transpositions could also prove interesting [START_REF] Müller | Proof-set search[END_REF][START_REF] Kishimoto | A solution to the GHI problem for depthfirst proof-number search[END_REF][START_REF] Saffidine | UCD: Upper Confidence bound for rooted Directed acyclic graphs[END_REF].

Fig. 1 .

 1 Fig. 1. Time needed to solve various opening positions in the game of y.

Table 2 .

 2 Selection policy for PP. C denotes the set of children.

	Node label Chosen child
	Max	arg max C PPN
	Min	arg min C PPN

Table 3 .

 3 Number of positions solved by each algorithm and number of positions on which each algorithm was performing best. PP, 4180 positions are solved fastest by MCTS, 30,719 positions need fewest iterations to be solved by PP, and 46,291 need fewest iterations by

		PP MCTS MCTS-E
	Positions solved 77,010 76,434 69,298
	Solved fastest	68,477 3,645	4,878
	Fewest iterations 22,621 35,444 18,942
	solved fastest by MCTS.

Table 4 .

 4 Number of node expansions needed to solve various sizes of domineering.

			5 × 6	6 × 6	7 × 6
	αβ	701,559 38,907,049 6,387,283,988
	PNT 1,002,277	>10 7	>10 7
	PN 2	17,236 >154,107	>511,568
	PP		101,244 5,525,608	>10 7
	PPT	27,766	528,032	4,294,785
	PP	2	3,634	24,190	145,757

Table 5 .

 5 Time (s) needed to solve various sizes of domineering.nogo is the misere version of the game of go. It was presented in the BIRS 2011 workshop on combinatorial game theory[START_REF] Chou | Revisiting Monte-Carlo tree search on a normal form game: Nogo[END_REF].3 The first player to capture has lost.We present results for the following algorithms: αβ, PNT[START_REF] Schijf | Proof-number search and transpositions[END_REF], PN 2 [5], PP, PPT and PP 2 . Again, the PNS algorithm could not find a single solution within 10 7 node expansion and is left out.

		5 × 6	6 × 6	7 × 6
	αβ 0.87	40.68	5,656
	PNT 5.92
	PN 2 78.7 >10,660 >153,000
	PP	0.24	20.1 >35.84
	PPT 0.17	5.33	55.13
	PP 2	0.22	15.5	320.3
	4.3 nogo

Table 6 .

 6 Number of node expansions needed to solve various sizes of nogo.

		4 × 4	18 × 1	20 × 1	22 × 1
	αβ 17,194,590 4,444,384 154,006,001 3,133,818,285
	PNT 3,575,076 2,015,179	>10 7	>10 7
	PN 2	77,010 > 22, 679 > 29, 098
	PP	>10 7 864,951 6,173,393	>10 7
	PPT 2,319,816	98,991	389,119	2,814,553
	PP 2		14,246

Table 7 .

 7 Time (s) needed to solve various sizes of nogo.

		4 × 4 18 × 1 20 × 1 22 × 1
	αβ	33.05	10.43	361.0	7,564
	PNT	436.6	144.2 > 809
	PN 2 27,519 > 3, 607 > 4, 583
	PP > 338.84	21.39	156.3 > 307.55
	PPT	396.36	9.46	46.3 446.58
	PP 2		109.7

Akihiro Kishimoto, personnal communication.

Some results can also be found on http://www.personeel.unimaas.nl/ uiterwijk/Domineering_results.html.

http://www.birs.ca/events/2011/5-day-workshops/11w5073