Grégoire Uhlrich
email: g.uhlrich@ipnl.in2p3.fr

Farvah Mahmoudi

Alexandre Arbey

A C++ framework automating theoretical calculations Beyond the Standard Model

Keywords: BSM, 1-loop calculation, Automated, Wilson coefficient

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.

Introduction

Calculations Beyond the Standard Model (BSM) have been a challenge for a long time, especially at the one-loop level. Transition amplitudes, differential cross-sections, and Wilson coefficients need to be evaluated in order to obtain the values for physical observables. Comparing them to the SM predictions and experimental results provides means to discriminate BSM models.

The one-loop order is often needed, as some relevant processes only appear at this order. This is the case of Flavour Changing Neutral Currents (FCNC) in flavour physics [START_REF] Buras | Weak Hamiltonian, CP Violation and Rare Decays[END_REF], that require one-loop values for Wilson coefficients. At the one-loop level, calculations cannot be done numerically as the very large number of terms together with renormalization processes require symbolic computations similar to what can be done by hand.

Currently, symbolic calculations at the one-loop level are done mainly using Mathematica [2], a commercial and closed computer algebra system. Several packages based on Mathematica implement high-energy physics calculations. FeynRules [3] computes Feynman rules from a BSM Lagrangian, from which FormCalc [START_REF] Hahn | Automatized one loop calculations in fourdimensions and D-dimensions[END_REF] can derive tree-level and one-loop quantities such as transition amplitudes or differential cross-sections. Packages such as Form-Flavor [START_REF] Evans | FormFlavor Manual[END_REF] or FlavorKit [START_REF] Staub | Exploring new models in all detail with SARAH[END_REF][START_REF] Porod | A Flavor Kit for BSM models[END_REF] calculate Wilson coefficients using the FormCalc machinery.

Many efforts have been employed to design independent computer algebra systems for high-energy physics. GiNaC [8] is a C++ library for symbolic computations written for that purpose. LanHEP [START_REF] Semenov | LanHEP: A Package for the automatic generation of Feynman rules in field theory. Version 3.0[END_REF], CompHEP [START_REF] Boos | CompHEP 4.4: Automatic computations from Lagrangians to events[END_REF] and CalcHEP [START_REF] Pukhov | CompHEP: A Package for evaluation of Feynman diagrams and integration over multiparticle phase space[END_REF] automate together calculations from the Lagrangian using their own symbolic computation framework. However, Mathematica-based packages are still the only ones able to provide one-loop calculations and Wilson coefficients for BSM.

We introduce MARTY, a modern and user-friendly solution to this issue. It provides for the first time a unique, free and open-source code fully written in modern C++ (2017 standard), implementing all the theoretical BSM machinery up to the one-loop order. Amplitudes, cross-sections and Wilson coefficients up to dimension 5 operators can be computed automatically. Support for dimension 6 operators will come in a future version. MARTY does not depend on any other framework, as it includes its own C++ Symbolic computation Library called CSL to manipulate mathematical expressions. It also comes with GRAFED, a desktop application displaying Feynman dia-grams on screen.

The reader can find the code, manuals, the documentation and more information on https://marty.in2p3.fr.

Code overview

MARTY is organised in several modules. They have been logically separated to ensure the code to be modular and general. The results are three modules: the physics core of MARTY, and two other modules which are both fully independent from each other and can be used as standalone, namely the computer algebra system of MARTY (CSL) and the application drawing Feynman diagrams (GRAFED). A schematic view of MARTY 's design is presented in figure 1.

MARTY (physics core)

This is the main part of the code, which contains all physics calculations, conventions, models, simplifications, etc. It uses its own CSL module as a mathematical backend to perform calculations specific to high-energy physics and the GRAFED module to draw Feynman diagrams. It also contains all group theory implementations, and model building features. Amplitudes, differential partonic cross-sections and Wilson coefficient may be calculated at tree-level or at the one-loop order. All calculations are automatic, symbolic, and can be performed in a very large variety of models as detailed in section 5.

CSL (computer algebra system)

This module does not know anything about physics, is logically separated from the physics part and can be used independently. It is a C++ Symbolic computation Library allowing us to handle mathematical expressions, tensors and simplifications needed to perform high-energy physics calculations. It is not as comprehensive as a standard computer algebra system like Mathematica because many features were not required for particle physics. It may however be extended in that direction if needed.

GRAFED (Feynman diagram generation and edition)

When doing calculations in particle physics, it is often convenient to visualize what the code is doing, and possibly include diagrams in publications. GRAFED was developed for this purpose and is also fully independent of the other modules of MARTY. It has three major features:

• An algorithm that finds an optimal way to place nodes in a 2D space to display Feynman diagrams. This algorithm is fully general (with no limit in the diagram size or number of loops) and automated. This allows one to quickly draw all diagrams for a particular process, without asking anything from the user and independently of the diagram topologies.

• A Graphical User Interface (GUI) that displays the generated diagrams. When asked, MARTY will run GRAFED with all the diagrams of a particular process. These diagrams appear then in the GUI, and may be exported (as png files or LATEX codes for the tikz-feynman package) directly to be included in a publication for example.

• The possibility to edit or create diagrams from scratch. Diagrams generated automatically by GRAFED are rather neat, but there is the possibility to edit, graphically, any aspects of the diagram (nodes, edges, labels, layout, etc) very easily. One can also create diagrams from scratch using GRAFED independently of MARTY.

MARTY design philosophy

The design of MARTY is guided by strong principles ensuring a final result corresponding to programming standards. First, the general principles unrelated to physics are:

• Independence. MARTY is written from scratch and is thus fully independent of any other framework. As such, there is no limit in what can be implemented in the code. Developers have a full control on any aspects of MARTY, to modify or extend its capabilities.

• User-friendliness. The code must be easy to use. The fact that it is written in C++ is a supplementary challenge in that purpose, but a modern knowledge of this language provides freedom for the userinterface. We think that this objective is fulfilled, since the normal usage of MARTY does not require any particular C++ knowledge and would be similar in many languages, including Mathematica.

• Modularity. MARTY is built as modular as possible. This means that unnecessary logical connections between different parts of the code are avoided. This is an important advantage for maintainability, since replacing or correcting a part of the code will become simpler.

• Readability. It is important for a code to be easily understandable by everyone, and in particular by a user willing to further develop the code. Strict coding conventions, clear naming for files / functions / variables and clever separation of different logical units make MARTY easy to understand considering its large size.

• Performance. C++ and python were the two main languages possible for MARTY as they are well-known in the high energy physics community. The choice of the language, C++, is related to performance reasons. A C++ code will run much faster in average than python for this type of code.

Concerning physics aspects, MARTY has been written with the following aims in mind:

• Generality. MARTY is designed to be as general as possible in the models it can handle, algebraic simplifications it can do, and calculations automated with it. A high level of generality has already been reached, and further developments will continue to focus on this aspect.

In particular there is no hard-coding because MARTY is expected to be extended even further in its future developments.

• Model independent calculations. In order to have an easy-to-use code for high energy physics, the computations have to be done in a model-independent way. The same code computing a given quantity should work for all models. Studying new models would hence imply only to write the Lagrangian or Feynman rules associated to it, and then using the same scripts to calculate the same quantities in the given model.

• One-loop level automated calculations. Calculations in BSM phenomenology often require at least one-loop level quantities. Many processes are trivial at tree-level but higher order corrections can be important from a phenomenological point of view when studying BSM models, as it is the case for instance for FCNC decays in flavour physics.

The one-loop level being significantly more difficult to calculate by hand than the tree-level, it is important to automate its calculation.

The efforts made to respect this philosophy will be useful in MARTY 's future developments as well. A code as general and independent as MARTY could benefit from a community effort to be maintained and developed and we think that the way it is written would allow for such a collaborative work.

Installation and usage

MARTY is available for download from its website: https://marty. in2p3.fr where one can find the manuals, documentation, publications and more information.

Installation

MARTY is open-source, GPL3 licenced and written in C++. By default, the installation script will use gcc-7, g++-7, gfortran-7 to compile MARTY and its dependencies. In case another version is more suitable, one can write

$ export MARTY_COMPILER_VERSION=n

with n the version number (if 0 is given the default compilers will be used instead but this is not recommended).

On Ubuntu1 / Debian since MARTY-1.0 (other Linux distributions may also be compatible but have not been tested) and on Mac since MARTY-1.1, the code can be built and installed using the following instructions: $ source setup.sh <installation-path> $ make $ make install

The installation path is /usr/local/ by default. Paths for header files, library files and binaries (the GRAFED application) must appear in environment variables in order to compile MARTY programs easily. This will be done by the setup.sh script for the current bash session, if MARTY is installed with the corresponding make target. This script must be executed using source to allow one to compile C++ programs using MARTY as presented in the following. Once the installation is performed, a marty_env.sh script is created in the root directory that contains the required bash instructions to set all environment variables. It can be executed at any time, again using source $ source marty_env.sh One can also copy its content in a .bashrc configuration file for example, to automate its execution for all bash sessions.

A procedure exists to use a MARTY or CSL program without installing it. A program.cpp file in marty/scripts can be compiled and executed typing from the main directory: The project layout specific to the marty/ directory together with the other modules (CSL and GRAFED) is presented in figure 2.

$ cd marty $ make program.x $ bin/program.x

Usage (scripting)

As a C++ framework, MARTY can be used in a C++ program after writing a few lines to include it, provided it is installed on the computer2 : #include <marty.h> using namespace mty; using namespace csl; using namespace std;

The three lines of using namespace are not necessary, but will allow us to omit prefixes mty:: (MARTY ) csl:: (CSL) and std:: (C++ standard library) in front of objects and functions. The main function containing the program can now be written:

int main() {
// The MARTY program return 0; }

Usage (compilation and execution)

In C++, source files need to be compiled before being run. MARTY uses the C++17 standard that appears in compiler options. To compile a source file main.cpp into an executable main.x the following commands are necessary 3 :

$ g++ -std=c++17 -c main.cpp -o main.o $ g++ -std=c++17 -o main.x main.o -lmarty $ ./main.x

Dependencies

MARTY has been written from scratch. Thus it has no dependency before the numerical evaluation of symbolic results. It contains in particular its own computer algebra system, CSL (C++ Symbolic computation Library), as a separate module. Therefore for the physics calculations from the Lagrangian to a final one-loop result simplified as much as possible, MARTY uses nothing but its own code and the C++ standard library.

For the numerical evaluation of results, there are two dependencies. The first, LoopTools [START_REF] Hahn | Automatized one loop calculations in fourdimensions and D-dimensions[END_REF], provides numerical values for scalar integrals arising at the one-loop level. Momentum integrals have always the same form, like the following 3-point function integral

I 3 ≡ d 4 q iπ 2 p µ 1 p ν 2 (q 2 -m 2 0 )((q + p 1 ) 2 -m 2 1 )((q + p 2 ) 2 -m 2 ) 2 .
(1)

3 g++ may be replaced by another C compiler, clang for example.

A way to treat this kind of integrals (including regularization) is to decompose the result in different possible Lorentz structures, each having a scalar factor in front, that is calculable with standard prescriptions [START_REF] Ellis | One-loop calculations in quantum field theory: From feynman diagrams to unitarity cuts[END_REF]. The decomposition for the 3-point function I 3 is the following

I 3 ≡ C 00 g µν + C 11 p µ 1 p ν 1 + C 12 (p µ 1 p ν 2 + p µ 2 p ν 1 ) + C 22 p µ 2 p ν 2 , (2)
with C ij scalar form factors depending on masses and squared momenta in the loop. The decomposition is done analytically by MARTY, and the evaluation uses LoopTools functions to determine the values of C ij .

The second dependency for the numerical evaluation is GSL [START_REF]Gnu scientific library[END_REF], a wellknown numerical library for C and C++. For complicated models with non trivial mixings (such as supersymmetric models), one has to diagonalize mass matrices to obtain the mass spectrum and mixings of the theory. For example from a non-diagonal squared-mass matrix M 2 of Φ, one calculates the eigenvector Φ that diagonalizes the matrix to

M 2 D Φ = U Φ , (3)
M 2 D = U † M 2 U. (4)
The symbolic result of MARTY is fully general (unless specified otherwise by the user) and uses generic symbols for all masses and mixings (matrices M 2 D and U in the example). For the numerical evaluation, input parameters must be given by the user. The diagonalization is then performed to get the spectrum and mixings, and to calculate the final results. The numerical diagonalization is performed using GSL.

Finally, GRAFED has a Graphical User Interface (GUI) that uses Qt [START_REF]The Qt Company, Qt open source model[END_REF]. Qt is a C++ framework allowing to build desktop applications fairly easily and is free and open-source with a GPL licence.

MARTY 's capabilities

This section presents in detail the calculations that can be performed with MARTY, the possible models, and the outputs that the code returns to the user.

Model building

First it is important to have a clear view of the BSM models that MARTY can handle. A model lies in a 4-dimensional Minkowski space-time and is defined by:

• A gauge group. The gauge group may be any combination of Semisimple Lie Groups4 (U (1), SU (N ), SO(N ),

Sp(N ), E 6 , E 7 , E 8 , F 4 , G 2 ). The unbroken Standard Model gauge SU (3) c × SU (2) L × U (1)
Y is an example of such a combination.

• A particle content. Each particle is an irreducible representation of the gauge group, i.e. an irreducible representation of each group composing the gauge. A particle may have a spin 0, 1/2 (Weyl, Dirac or Majorana), or 1. All gauge couplings are introduced automatically by MARTY without any help from the user.

• Additional couplings. The user can add any interaction term in the Lagrangian. MARTY simply checks that combining the unbroken gauge representations of the interacting particles gives indeed a trivial representation 5 .

There are two ways to build a model in MARTY. The first one is the most straightforward way but also the most complicated one. It consists in giving explicitly the full Lagrangian to MARTY. Few terms in general are provided by unbroken gauge couplings, in particular when one studies a phenomenological model extending the SM. In the SM, there are about 100 terms to write by hand coming from the symmetry breaking. In the Minimal Supersymmetric extension of the Standard Model (MSSM), several thousands. It is possible to do it but one has to be very careful on every convention, sign and factor in front of each term. A small error can lead to wrong results due to interference between different diagrams for a given process.

The second option is to define a high energy Lagrangian with all symmetries preserved, and give MARTY prescriptions to break it. The initial Lagrangian has much less interaction terms and is simpler to write. Based on correct prescriptions (gauge, flavour symmetry breaking, replacements, renaming, etc) MARTY will basically re-derive the final Lagrangian for the user. This solution will not necessarily be the easiest one depending on the model but is certainly a practical option. It is in particular the way chosen to build the MSSM in MARTY.

In the following we present a sample code building a SU (2) L gauge with one quark in the doublet representation

Q L = u L d L , (5)
broken by MARTY with a single instruction. We also ask the code to rename the broken fermions Q 1 and Q 2 to u and d which corresponds to standard conventions. With 3 lines of breaking prescriptions MARTY derives the 17 interaction terms (including vector-ghost interactions 6 ) between the final 8 particles in the model. For a more evolved model one needs more instructions to specify every conventions, but this method is still very practical and more efficient than giving the full Lagrangian by hand. For a gauge symmetry breaking, explicit expressions of gauge generators (T A ij , f ABC ) must be known and MARTY will not necessarily know them. For SU (2) and SU (3) SM gauge terms, the whole procedure is automated, but for other groups the user may have to define expressions for generators.

MARTY also contains built-in models that can be used directly for calculations: Scalar φ 3 theory, Scalar QED, QED, QCD, Electroweak model, Standard Model, 2 Higgs Doublet Models and Minimal Supersymmetric Standard Models (unconstrained and phenomenological).

Amplitudes

Transition amplitudes from an initial state to a final state noted iM(i → f ) are the basic quantities that MARTY is able to calculate. It uses the Lagrangian exponentiation as well as the Wick's theorem [START_REF] Wick | The evaluation of the collision matrix[END_REF] to find all possible diagrams and derive their corresponding expressions. This step is fully general and has no limit in the diagram complexity or in the number of external legs. Amplitudes are then used to calculate differential cross-sections or to derive Wilson coefficients.

Once an analytical expression for a given diagram has been found, it needs to be simplified in several ways in order to obtain a numerical evaluation of the result. Simplification steps done by MARTY are the following: Dirac algebra simplification. This includes calculation of traces in the Dirac space and simplifications in γ-matrix products, for particles of spin 1/2 [START_REF] Schwartz | Quantum Field Theory and the Standard Model[END_REF].

Group algebra simplification. Similarly to γ-matrices, algebra generators have to be simplified into amplitudes. Projection operators are used [START_REF] Cvitanovic | Group theory for Feynman diagrams in non-Abelian gauge theories[END_REF] and traces are calculated in all semi-simple groups [START_REF] Van Ritbergen | Group theory factors for Feynman diagrams[END_REF]. The remaining colour structures that cannot be simplified are stored and factored from the rest of the amplitude, in dedicated abbreviations. For standard gauge groups and in particular for fundamental representations all possible terms will be simplified automatically 7 .

Minkowski Index contraction. Minkowski indices are expanded and contracted as much as possible in D-dimensions to perform Dimensional Regularization (DREG) at the one-loop order. As in D-dimension g µ µ = D, one has to expand the whole diagram to gather all factors of D.

Reduction of one-loop momentum integrals.

A momentum integral at oneloop can be decomposed on the basis of scalar form factors [START_REF] Ellis | One-loop calculations in quantum field theory: From feynman diagrams to unitarity cuts[END_REF]. These form factors depend on masses and momenta and can be provided by LoopTools [START_REF] Hahn | Automatized one loop calculations in fourdimensions and D-dimensions[END_REF] up to the rank 4 5-point functions, i.e. loops with five external legs and four momenta in the numerator. This is the actual limit for fully-simplified one-loop quantities in MARTY.

Dimensional Regularization. The form factors coming from one-loop integrals can have a divergent part that is regularized by taking the dimension D = 4 -2 . In this case, integrals take the form

I ≈ a + b + O( ). (6)
Factors of D coming from Minkowski index contractions must then be kept to determine the local terms they generate when they are multiplied by a divergent integral [START_REF] Denner | Techniques for the calculation of electroweak radiative corrections at the one-loop level and results for W -physics at[END_REF][START_REF] Sulyok | A closed expression for the uv-divergent parts of one-loop tensor integrals in dimensional regularization[END_REF]. For the scalar one-point function for example, we get the finite part

Finite(DA 0 (m 2 )) = Finite((4-2 )A 0 (m 2 )) = -2m 2 +4•Finite(A 0 (m 2 )). (7)
Equations of motion. For spin 1 particles, the equation of motion is simply

µ (p)p µ = 0, (8)
where (p) is the polarization 4-vector of the boson. For spin 1/2 particles, the Dirac equation is applied. It reads

/ pu(p) = mu(p) for particles, (9)
v(p) / p = -mv(p) for anti-particles.

Factorization. Results are partially factored to compactify at most the final expressions. In particular, factorization by masses and momenta are performed.

Abbreviation. Abbreviations are introduced automatically by MARTY to lighten expressions and gain in execution time. All abbreviations used can be displayed by typing DisplayAbbreviations().

Using the toy model presented in section 5.1, the calculation of an amplitude in MARTY is very simple. One has to give the order (mty::Order::TreeLevel here), the model, and the field insertions. Let us consider the process uū → d d at tree-level, the following instruction is needed to run the computation:

auto res = ComputeAmplitude( Order::TreeLevel, model, {Incoming("u"), Incoming(AntiPart("u")), Outgoing("d"), Outgoing(AntiPart("d"))} );

The code is in C++, therefore variable types must be declared. However since C++11, the compiler does not need the user to specify the type returned by a function anymore. By typing auto the compiler is told to find the exact type itself. Then, res.expressions contains the different terms of the amplitude and res.diagrams the Feynman diagrams. To display expressions in standard output, show the diagrams in GRAFED and display the abbreviations introduced by MARTY the following three lines are necessary respectively: Display(res); Show(res); DisplayAbbreviations(); MARTY 's output is presented in the following, displaying diagrams in GRAFED as the screenshot in figure 3 shows. Only the index IDs have been changed to lighten the output: 0 : 1/8*Ab_0001*EXT_{k,\%eps_1,i,\%del_1,l,\%del_2,j,\%eps_2} *gamma_{+\%sigma,\%del_2,\%eps_1}*gamma_{\%sigma,\%del_1,\%eps_2} 1 : 1/8*Ab_0001*EXT_{k,\%eps_1,i,\%del_1,l,\%del_2,j,\%eps_2} *gamma_{+\%sigma,\%del_2,\%eps_1}*gamma_{\%sigma,\%del_1,\%eps_2} 2 : -1/8*Ab_0002*EXT_{k,\%eps_1,i,\%del_1,l,\%del_2,5,\%eps_2} *gamma_{+\%sigma,\%del_1,\%eps_1}*gamma_{\%sigma,\%del_2,\%eps_2}

Total : 3 particle amplitudes.

The abbreviations read: The resulting amplitude can be read off form MARTY 's output, keeping a separate term for each diagram:

iM = ig 2 L 8 1 s 13 d(p 4 )γ µ u(p 1 )ū(p 2 )γ µ d(p 3 ) + 1 s 13 d(p 4 )γ µ u(p 1 )ū(p 2 )γ µ d(p 3 ) - 1 s 12 ū(p 2 )γ µ u(p 1 ) d(p 4 )γ µ d(p 3 ) , (11)
with s ij ≡ p i • p j .

Cross-sections

Cross-sections are the main observables used in collider physics. They are directly proportional to the number of events observed in the various detectors. MARTY does not compute directly the cross-sections but calculates the complicated theoretical part namely the squared amplitudes. For incoming particles {I} of spins {j I } and outgoing particles {O} of spins {j O } the squared amplitude is (as a function of the amplitude iM that depends on the particle spins) 1

I d I {j I },{j O } |M| 2 , ( 12
)
with d I the spin dimension of the incoming particle I taking into account massless effects for spin 1 particles. This quantity is averaged (summed) over the spin dimensions of incoming (outgoing) particles. Calculating the squared amplitudes implies the calculation of traces in Dirac and colour spaces (group algebra) that MARTY computes automatically. The result is a scalar depending on momenta and masses of particles in the process. The differential cross-section has always the same form for a given process of amplitude iM dσ ≡ K(p i , m i ) • 1

I d I {j I },{j O } |M| 2 dΠ LIP S , (13)
with K(p i , m i ) a factor coming from kinematics, and dΠ LIP S the Lorentz Invariant Phase Space. Once the amplitude squared has been calculated and simplified, no more computer algebra system is needed to pursue the calculation. This is the quantity that MARTY can compute automatically.

Considering the toy model of section 5.1, calculating the squared amplitude is very simple. The user must first calculate the amplitude, and simply square it with MARTY. The average over incoming spins is done by MARTY, i.e. the returned quantity corresponds to equation 12. After calculating this quantity with MARTY, the user has again one single line to write:

Expr squared_ampl = model.computeSquaredAmplitude(res); cout << "<|M|^2>␣=␣" << squared_ampl << endl; DisplayAbbreviations(); Expr is the main variable type of CSL, internal representation of a symbolic mathematical expression. The output in terminal is presented in the following.

<|M|^2> = 1/4*s_14*s_23*( 1/2*Ab_0001^(*)*Ab_0002 + Ab_0001*Ab_0001^(*) + 1/2*Ab_0001*Ab_0002^(*) + 1/4*Ab_0002*Ab_0002^(*)) Ab_0001 = i*g_L^2/s_13 Ab_0002 = i*g_L^2/s_12
One can see that abbreviations have been introduced by MARTY. They can be expanded and the result can be further factored by CSL typing // Evaluate abbreviations Evaluate(squared_ampl, eval::abbreviation); // Factor the whole expression DeepFactor(squared_ampl); cout << "<|M|^2>␣=␣" << squared_ampl << endl;

In this way, one can obtain a compact result:

<|M|^2> = 1/16*g_L^4*(s_12^(-2) + 4*s_13^(-2) + 4/(s_12*s_13)) *s_14*s_23
As one can see above, MARTY 's outputs contain scalar products of external momenta, namely s ij ≡ p i •p j . MARTY does not perform any kinematics for now, i.e. stops the simplification as in the output shown. This could be easily implemented in the future, for example introducing Mandelstam variables. However as it does not represent an important analytical challenge by hand and that it has no impact on the following numerical evaluation (see section 5.5), this part is for now left to the user.

Wilson coefficients

Wilson coefficients are the coefficients in front of particular operator structures in an amplitude [START_REF] Buras | Weak Hamiltonian, CP Violation and Rare Decays[END_REF]. For the b → sγ process that will be detailed in section 6, the amplitude may be decomposed on a two operator basis, each one with a scalar coefficient in front. Naming q the photon momentum and its polarization vector, one obtains

iM(b → sγ) = -4G F √ 2 e 16π 2 V tb V * ts m b C 7 Ô7 + C 7 Ô 7 , (14)
with

Ô( ) 7 ≡ sγ| Ô( ) 7 |b = sσ µν P R(L) bF µν , (15)
σ µν = i 2 [γ µ , γ ν ] , (16)
F µν = iq [µ ν] = i 2 (q µ ν -q ν µ ) . (17)
The global factor

-4G F √ 2 e 16π 2 V tb V *
ts m b is defined by convention. This procedure to decompose amplitudes in Wilson coefficients and operator matrix elements is used in particular in flavour physics. As quarks appear only in bound states, the partonic amplitude is not the full story. One has to take into account long-distance effects that cannot be calculated pertubatively. The b → sγ transition may correspond for example to a hadronic process B0 → K0 γ. These long-distance effects are model-independent and arise only in the operator matrix element between final and initial states F | Ô|I . The BSM dependence is then contained in the Wilson coefficient that can be calculated perturbatively by MARTY.

For now, operators of dimension 6 with 4 fermions cannot be directly given by MARTY at one-loop as some simplifications are needed that are not yet implemented. The step missing is a double application of Fiertz identities, to simplify all possible momenta in fermion currents. A Wilson coefficient for such an operator could still be read off in the results but would ask the user to do some algebra by hand, to determine which part of the amplitude contributes to the coefficient.

A concrete example of Wilson coefficient calculation is presented in section 6 which is devoted to the calculation of C 7 in the MSSM.

Library generation

Results of MARTY can in general not be used directly. Depending on the process, it may be a very big and complicated analytical expression. What the user may need are numbers, i.e. numerical evaluations of the analytical int main() { cout << "XSec␣=␣" << squared_ampl(0.1, 100, 60, 40, 40) << endl; return 0; } A library may contain as many functions as wanted. This procedure is fully general and is automated. Note that if the library needs additional include or library paths (in particular if CSL and MARTY are not installed in standard locations), it is possible to specify them with: myLib.addIPath("/home/.local/include"); myLib.addLPath("/home/.local/lib"); 5.6. Feynman diagrams GRAFED is the part of MARTY generating and rendering Feynman diagrams. It is used to create automatically diagrams when calculating a process with the Show(res) command as we discussed in section 5.2. It can also be used to edit or create diagrams from scratch. Many aspects of the diagrams can be chosen by the user in an intuitive way. A screenshot of GRAFED is shown in figure 4.

GRAFED will in the future be released as standalone. All diagrams in this publication are generated automatically or created with GRAFED.

6. Calculation of δ LO C χ, t 7 (M W ) in the pMSSM An example of MARTY 's capabilities is presented in this section, namely the calculation of the MSSM contribution to the Wilson coefficient C 7 . This coefficient describes the b → sγ transition and is non zero only at the oneloop level.

The pMSSM

We consider here the phenomenological Minimal Supersymmetric Standard Model (pMSSM) which is a generic CP conserving MSSM framework that has 19 parameters more than in the SM as opposed to the full MSSM which has 105 extra parameters. The pMSSM has been chosen for validation because of its complexity and generality. Obtaining a correct result in this model demonstrates MARTY 's capabilities for model building and symbolic calculations.

The Wilson coefficient C 7

We calculate the Leading Order (LO) value of the Wilson coefficient C 7 , associated with the operator in equation 15. The process is showed in figure 5. It is a FCNC process with a photon changing the quark flavour b → s which is forbidden at tree-level in the SM and the pMSSM, and the LO is thus at the one-loop level. Strong experimental constraints exist for FCNCs and their calculations for BSM models represent an important task for phenomenology. We consider in this example one of the supersymmetric contributions i.e. diagrams with the top squarks and charginos shown in figure 6. We perform the calculation on-shell in the Feynman-'t Hooft gauge 8 . The reversal of the fermion-flow in the diagram is due to fermion-number violating interactions between charginos and SM fermions. This is mostly related to the definition of charginos and may be treated following prescriptions of [START_REF] Denner | Compact Feynman rules for Majorana fermions[END_REF]. At the end of the calculation, the fermion flow is regular but may get a sign due to charge conjugation matrix C which appears. This sign is important to determine exactly because of interferences between the diagrams.

We vary two pMSSM parameters, µ (the Higgsino parameter) and M 2 (the Wino mass). More details on MSSM parameters are given in [START_REF] Martin | A Supersymmetry primer[END_REF]. Contributions to C 7 come from chargino and stop loops and depend on µ and M 2 in particular. The chargino mass matrix reads

M χ = 0 X T X 0 , (18)
with

X = M 2 √ 2 sin βM W √ 2 cos βM W µ , (19)
tan β being the angle between the two Higgs doublets' Vacuum Expectation Values (VEVs). with 10 identical digits in average, for both spectra. One explanation may be that we used quadruple precision (128 bits) floating point variables for FormCalc and MARTY 's outputs, and only double precision (64 bits) for SuperIso output and the analytical formula. The 4-digits precision is however completely satisfactory considering the uncertainty coming form higher orders in perturbation theory. This example completes the presentation of what MARTY can calculate. MARTY will generate for any process, in any model, libraries evaluating theoretical quantities and give the user a spectrum generator at the same time. In the case of supersymmetry, spectrum generators already exist with higher-order terms but in a general BSM models one needs this generic treelevel spectrum generator. More information and example can be found on the website https://marty.in2p3.fr.

Performance

We measure the performance of a computer program with two main indicators. The execution speed and the quantity of memory (RAM) the program needs for running. For BSM symbolic calculations at one-loop, it is not possible to give a standard execution speed nor the quantity of memory as it depends on the model and the process to calculate. The amount of memory taken by MARTY is typically very small. It is very rare to reach 1 GB, and is often under 100 MB. Indicative values of execution times are shown in table 2, measured on various processes, always running on a single CPU.

External legs Tree-level One-loop 2 ≤ 10 -1 10 -2 3 ≤ 10 -1 10 0 /10 It can be seen in table 2 that the calculation complexity depends strongly on the number of legs at one-loop. The more legs there is for the loop, the more terms appear in the amplitude. This number of terms grows very fast with the number of legs connected to the loop and explains the results shown here.

For squared amplitudes there is no simple rule to determine the execution time but it is in general several orders of magnitude more than the simple amplitude calculation as squaring the amplitude also squares the number of terms to simplify. Improving performance for this calculation is an important development for the next release of MARTY.

Future developments

MARTY has fulfilled most of the planned requirements, but further developments are ongoing which are listed in the following.

• Wilson coefficients for 4-fermion operators. For now, MARTY can compute amplitudes for 4-fermions processes, but cannot give automatically the corresponding Wilson coefficients because of a missing simplification step. This step is the double application of Fiertz identities to simplify all momenta in quark currents. Once this simplification will be implemented, 4-quarks operators will become available in MARTY.

• More group theory simplifications. All simplifications with algebra generators are not implemented in MARTY. Some are missing because it is very difficult to automate these identities for all semisimple groups and all representations. These missing simplifications concern mostly non-fundamental representations, exceptional algebras and squared amplitudes for pure gluonic amplitudes. Further developments will focus on this issue but the user can also define easily the missing properties.

• Automated NLO corrections. With MARTY one can calculate all one-loop quantities needed to renormalize a BSM model. However this procedure is not automated and will surely be a point of attention in the future.

• Operator mixing for Wilson coefficients. Renormalization comes with operator mixing for Wilson coefficients. This task is more challenging but there is currently no code able to fully automate this procedure for general BSM. Therefore having a code able to do this task would be very useful for flavour physics.

• Interfaces with other codes. This publication presents the first version of MARTY, that is at the moment not interfaced with other codes. This is the most important improvement that is planned for next developments. Universal Feynman Rules Output (UFO [START_REF] Degrande | UFO -The Universal FeynRules Output[END_REF]) will link MARTY directly to event generators, and work is in progress to have direct interface with SuperIso [START_REF] Mahmoudi | SuperIso: A Program for calculating the isospin asymmetry of B → K * γ in the MSSM[END_REF][START_REF] Mahmoudi | SuperIso v2.3: A Program for calculating flavor physics observables in Supersymmetry[END_REF][START_REF] Mahmoudi | SuperIso v3.0, flavor physics observables calculations: Extension to NMSSM[END_REF] for flavour physics, and SuperIso Relic [START_REF] Arbey | SuperIso Relic: A Program for calculating relic density and flavor physics observables in Supersymmetry[END_REF][START_REF] Arbey | SuperIso Relic v3.0: A program for calculating relic density and flavour physics observables: Extension to NMSSM[END_REF][START_REF] Arbey | SuperIso Relic v4: A program for calculating dark matter and flavour physics observables in Supersymmetry[END_REF] for dark matter phenomenology.

Conclusion and Outlook

We presented MARTY, a new C++ framework automating theoretical calculations symbolically for BSM physics. The degree of generality reached by MARTY has never been achieved before. It has its own computer algebra system (CSL) and automates all theoretical calculations directly from the Lagrangian. Feynman rules, Feynman diagrams, amplitudes, cross-sections, and Wilson coefficients can be obtained in a very large variety of BSM models up to the one-loop level. A full NLO treatment will also be implemented in the near future, treating renormalization of fields, masses, couplings, and Wilson coefficients (including operator mixings).

A proof of its capabilities has been demonstrated through a tree-level cross-section calculation, and a one-loop Wilson coefficient in the pMSSM. The results are at first symbolic mathematical expressions, but numerical C++ libraries are built automatically by MARTY allowing us to explore in full generality the parameter space of the model for some user-defined quantities. A spectrum generator specific to the user's model is also created automatically by MARTY when needed. Most of popular BSM models can be built in MARTY. The MSSM, extended gauge models, and vector-like quarks are examples of possible BSM implementations.

MARTY can already be very useful for BSM phenomenology in this current version. The particular advantage of MARTY is to be written as a unique code, not depending on any framework. Within MARTY, every aspect of model building and high-energy physics calculation are under control, in the same program and in the same language. This is a unique opportunity for future collaborations to take this code even further, extending it to new models, other simplification methods, or even different types of calculations.

Figure 1 :

 1 Figure 1: Schematic design of MARTY, including the three sub-modules: MARTY (physics part), CSL (symbolic manipulation) and GRAFED (Feynman diagram generation).

Figure 2 :

 2 Figure 2: Layout of the MARTY project. The main sub-directories are shown with their content.

  EXT = d_{k,del}(p_4)*d_{i,gam}(p_3)^(*) *u_{l,beta}(p_2)^(*)*u_{j,alpha}(p_1) Ab_0001 = i*g_L^2/s_13 Ab_0002 = i*g_L^2/s_12

Figure 3 :

 3 Figure 3: Screenshot of GRAFED displaying the three Feynman diagrams for the process uū → d d at the tree-level in the toy model defined in section 5.1. The tree diagrams correspond to the exchange of each broken vector from the SU (2) L initial symmetry, noted A Li by default in MARTY.

Figure 5 :

 5 Figure 5: b → sγ process represented in a model-independent way. The transition amplitude is the sum of all diagrams that can fill correctly the hatched disk. This diagram has been built using GRAFED.

Figure 6 :

 6 Figure 6: Two types of contribution for C 7 in the pMSSM, with stops t and charginos χ. These diagrams have been built using GRAFED.

Figure 7 :

 7 Figure 7: Results for C 7 (chargino and stop contributions) in the pMSSM, from MARTY one the left and from the analytical formula [23] on the right, for the spectrum generated by MARTY at tree-level. The results match to four digits in average.

Figure 8 :

 8 Figure8: Results for C 7 (chargino and stop contributions) in the pMSSM, from MARTY one the left and from the output of SuperIso[START_REF] Mahmoudi | SuperIso v2.3: A Program for calculating flavor physics observables in Supersymmetry[END_REF] on the right, for the spectrum generated by SOFTSUSY[START_REF] Allanach | SOFTSUSY: a program for calculating supersymmetric spectra[END_REF][START_REF] Allanach | The Inclusion of Two-Loop SUSYQCD Corrections to Gluino and Squark Pole Masses in the Minimal and Next-to-Minimal Supersymmetric Standard Model: SOFTSUSY3.7[END_REF] with two-loop corrections. The results match to four digits in average.

Table 2 :

 2 Typical execution time (order of magnitude) of MARTY in second per 100 Feynman diagrams during the calculation of an amplitude, for different numbers of external legs.

	1

On Ubuntu 16.04 or prior versions, MARTY can be built properly only since version 1.2.

For this to work properly, one needs either to install MARTY in the standard location (by default, /usr) or make sure that the paths to access include files in C++ (CPATH environment variable on Ubuntu), libraries (LIBRARY_PATH and LD_LIBRARY_PATH environment variable on Ubuntu) and binaries (PATH) contain the installation path of MARTY.

Strictly speaking, groups that have a semi-simple Lie algebra.

If it does not, this is the sign of an obvious gauge violating term. If it does, the term may still violate gauge symmetry but it is more difficult to test automatically.

Gauge fixing terms for ghosts are not taken into account here.

Group simplifications are easy to implement in MARTY in case there are missing ones. The complicated part is to determine theoretically the simplification rules to apply.

Other gauges can be used such as the unitary or Lorentz gauges.

project/5fb7c30c3175f910d2392059 CPC Library link to program files: (to be added by Technical Editor) Code Ocean capsule: (to be added by Technical Editor) Licensing provisions: GPLv3

https: // marty. in2p3. fr https://fr.overleaf.com/

results, for a given set of values of the model parameters. Let us consider the cross-section of section 5.3. The result is rather simple, but the principle would be exactly the same for a more complicated expression. The way it works in MARTY is also rather simple and is contained in a few lines (giving a library name and a path to create it): mty::Library myLib("uubar_to_ddbar", "."); myLib.addFunction("squared_ampl", squared_ampl); myLib.build();

A mty::Library is an abstract object that takes symbolic expressions as functions (here the cross-section), creates and compiles a C++ library allowing to evaluate them numerically. The function generated by MARTY is: The function takes as arguments all symbols (possibly complex here) that did not contain any value at the time the library was generated. The library is compiled automatically and can be used as demonstrated in the following.

#include "uubar_to_ddbar.h" using namespace std; using namespace uubar_to_ddbar;

The stop squared mass matrix reads

m Q 3 , m u 3 are soft supersymmetry breaking parameters, A t is a trilinear coupling, y t the top Yukawa and finally

The exact numerical values of pMSSM parameters used to evaluate C 7 are presented in table 1.

[-1000, 1000] GeV The results are shown in figures 7 and 8. MARTY 's output is compared with the analytical formula given in [START_REF] Ciuchini | Next-to-leading QCD corrections to B → X s γ in supersymmetry[END_REF] and with SuperIso [START_REF] Mahmoudi | SuperIso: A Program for calculating the isospin asymmetry of B → K * γ in the MSSM[END_REF][START_REF] Mahmoudi | SuperIso v2.3: A Program for calculating flavor physics observables in Supersymmetry[END_REF][START_REF] Mahmoudi | SuperIso v3.0, flavor physics observables calculations: Extension to NMSSM[END_REF]. Numerical evaluations have been done for two different spectra. This first one (figure 7) is a tree-level spectrum computed by MARTY using GSL [START_REF]Gnu scientific library[END_REF], and the result is compared with the analytical formula. The second spectrum (figure 8) is calculated by SOFTSUSY [START_REF] Allanach | SOFTSUSY: a program for calculating supersymmetric spectra[END_REF][START_REF] Allanach | The Inclusion of Two-Loop SUSYQCD Corrections to Gluino and Squark Pole Masses in the Minimal and Next-to-Minimal Supersymmetric Standard Model: SOFTSUSY3.7[END_REF] with two-loop order corrections which are known to be important for the charginos [START_REF] Martin | A Supersymmetry primer[END_REF]. For this spectrum, we compare MARTY with the output of SuperIso.

As can be seen for all the results there exists an excellent agreement between the analytical formula, SuperIso, and MARTY. The agreement is up to 4 digits. In addition, we also tested the results given by FormCalc [START_REF] Hahn | Automatized one loop calculations in fourdimensions and D-dimensions[END_REF] for the same process, and there is in this case a perfect agreement with MARTY


Vinh-Thong Ta
email: vinhthong.ta@unicaen.fr

Abderrahim Elmoataz
email: abderrahim.elmoataz-billah@unicaen.fr

Olivier Lézoray
email: olivier.lezoray@unicaen.fr

Adaptation of Eikonal Equation over Weighted Graph

In this paper, an adaptation of the eikonal equation is proposed by considering the latter on weighted graphs of arbitrary structure. This novel approach is based on a family of discrete morphological local and nonlocal gradients expressed by partial difference equations (PdEs). Our formulation of the eikonal equation on weighted graphs generalizes local and nonlocal configurations in the context of image processing and extends this equation for the processing of any unorganized high dimensional discrete data that can be represented by a graph. Our approach leads to a unified formulation for image segmentation and high dimensional irregular data processing.

Introduction

Solutions of the nonlinear eikonal equation have found numerous applications. One can quote for instance, geometric optics, image analysis or computer vision including shape from shading [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF][START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF], median axis or skeleton extraction [START_REF] Siddiqi | The hamilton-jacobi skeleton[END_REF], topographic segmentation (watershed) [START_REF] Maragos | Curve evolution, differential morphology and distance transforms as applied to multiscale and eikonal problems[END_REF] or geodesic distance computation on discrete and parametric surfaces [START_REF] Bronstein | Weighted distance maps computation on parametric three-dimensional manifolds[END_REF][START_REF] Sethian | Ordered upwind methods for static hamiton-jacabi equations: Theory and algorithms[END_REF][START_REF] Abgrall | Numerical discretization of the first-order hamilton-jacobi equations on triangular meshes[END_REF][START_REF] Shu | High order WENO schemes for hamiltonjacobi equations on triangular meshes[END_REF][START_REF] Mémoli | Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces[END_REF]. The latter works consider both structured and unstructured meshes on cartesian or non-cartesian domains.

The eikonal equation is a special case of the following general continuous Hamilton-Jabobi equation:

H(x, f, ∇f ) = 0 x∈Ω ⊂ IR n f (x) = φ(x) x∈Γ ⊂ Ω , (1)
where φ in the boundary condition is a positive speed function defined on Ω and f (x) is the traveling time or distance from source Γ . Then, the eikonal equation can be expressed by using the following Hamiltonian:

H(x, f, ∇f ) = ∇f (x) -P (x) , (2)
where P (x) is a given potential function. Solution of (1) represents the shortest distance from x to the zero distance curve given by Γ (where φ(x)=0). Solutions of (2) are usually based on a discretization of the Hamiltonian where the approximation of the derivatives is performed by the Godunov [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] or the Lax-Friedrich [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] schemes. Then, many numerical methods have been proposed and investigated to solve the nonlinear system described by [START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF]. For instance, one can quote the following schemes. (i) An iterative scheme [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] relying on fixed point methods that solves a quadratic equation was proposed. (ii) The fast sweeping methods [START_REF] Zhao | Fast sweeping method for eikonal equations[END_REF] that use Gauss-Seidel type of iterations to update the distance function field. The key point of fast sweeping is to update the points in a certain order. (iii) Tsitsiklis [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] was the first to develop a Dijkstra like method and proposed an optimal algorithm for solving the eikonal equation. Based on this idea, [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] produced the fast marching methods.

Another approach to solve (2) is to consider a time dependent version of the equation and to evolve it to the steady state. Then, (2) can be rewritten as

     ∂f (x, t)/∂t = -∇f (x) + P (x) x∈Ω ⊂ IR n f (x, t) = φ(x) x∈Γ ⊂ IR n f (x, 0) = φ 0 (x) x∈Ω . (3)
This paper only considers the discrete analogue of the time dependent formulation of the eikonal equation but, in future works, the stationary case (time independent) will be also considered.

Contributions.

In this work, we propose an adaptation of (3) over weighted graphs of the arbitrary structure. The goal here is to provide a simple and common formulation that solves the eikonal equation for any discrete data that can be represented by a weighted graph such as images or high dimensional data defined on irregular domains. This alternative formulation for solving the eikonal equation is based on partial difference equations (PdEs) and discrete gradients over weighted graphs. Our formulation has several advantages. Any discrete domain that can be described by a graph can be considered without any spatial discretization. In the context of image processing, local and nonlocal configurations are directly enabled within a same formulation. Finally, the aim of this paper is not to solve a particular application with the eikonal equation but to show the potentialities of our proposition to address image segmentation, data clustering or distance computation.

Paper Organization. The paper is organized as follows. Section 2 recalls basics, definitions and operators on weighted graphs. Section 3 introduces our formulation for solving the eikonal equation. Section 4 shows the potentialities of our proposition for the segmentation of images and unorganized data processing. Finally, last Section concludes.

Discrete Derivatives on Weighted Graphs

This Section recalls basics, definitions, operators and processes on weighted graphs.

Definitions and Weighted Graphs Construction

Notations and Definitions. We consider the general situation where any discrete domain can be viewed as a weighted graph. Let G=(V, E, w) be a weighted graph composed of two finites sets: vertices V and weighted edges E⊆V ×V . An edge (u, v)∈E connects two adjacent (neighbor) vertices u and v . The neighborhood of a vertex u is noted N (u)={v∈V \{u} : (u, v)∈E}. The weight ω uv of an edge (u, v) can be defined with a function w:V ×V →IR + such that w(u, v)=ω uv if (u, v)∈E and w(u, v)=0 otherwise. Graphs are assumed to be simple, connected and undirected implying that function w is symmetric. Let f :V →IR be a discrete real-valued function that assigns a real value f (u) to each vertex u∈V . We denote by H(V ) the Hilbert space of such functions defined on V .

Weighted Graphs Construction. Any discrete domain can be represented by a weighted graph where functions of H(V ) represents the data to process. In the general case, an unorganized set of points V ⊂IR n can be seen as a function f 0 :V ⊂IR n →IR n . Then, constructing a graph from this data consists in defining the set of edges E by modeling the neighborhood. It is based on a similarity relationship between data with a pairwise distance measure µ:V ×V →IR + . There exists several methods to transform a set of vertices V into a neighborhood (similarity) graph (see [15] for a survey on proximity and neighborhood graphs). In this paper, we focus on two particular graphs: the τ -neighborhood graphs and a modified version of k-nearest neighbors graphs. The k nearest neighbors graph, noted k-NNG is a weighted graph where each vertex u∈V is connected to its k nearest neighbors which have the smallest distance measure towards u according to function µ. Since this graph is directed, a modified version of this graph is used to make it undirected. The τ -neighborhood graph, noted G τ is a weighted graph where the τ -neighborhood N τ for a given vertex u∈V is defined as N τ (u)={v∈V \{u} : µ(u, v)≤τ } with τ >0 a threshold parameter.

2D images can be viewed as functions f 0 :V ⊂Z Z 2 →IR n . In this case, the associated distance µ for construct the neighborhood graph is usually the city block or the Chebychev distances computed with the spatial coordinates of each vertex representing an image pixel. With these distances and the τ -neighborhood graphs, one recovers the two usual graphs used in image processing, the 4adjacency grid graph (denoted G 0 with the city block distance) and the 8adjacency grid graph (denoted G 1 with the Chebychev distance) with τ ≤1. Another useful graph structure in image processing is the region adjacency graph (RAG) where vertices correspond to image regions, and the set of edges is obtained by considering an adjacency distance. With the τ -neighborhood (τ =1), the RAG is the Delaunay graph of an image partition.

Weights Computation. Similarities between data can be incorporated within edges' weights according to a measure of similarity g:E→IR + with w(u, v)=g(u, v) for (u, v)∈E. Then, the distance computation between data is performed by comparing their features that generally depend on a given initial function f 0 ∈H(V ). To this aim, each vertex u∈V is assigned with a feature vector F (f 0 , u)∈IR m . With F , the following weight functions can be considered. For a given edge (u, v)∈E and a distance measure ρ:V ×V →IR + associated to F , we can have

g 0 (u, v) = 1 (constant weight case) , g 1 (u, v) = (ρ(F (f 0 , u), F (f 0 , v)) + ǫ) -1 with ǫ>0, ǫ→0 , g 2 (u, v) = exp(-ρ(F (f 0 , u), F (f 0 , v)) 2 /σ 2 ) with σ>0 ,
where σ controls the similarity and ρ is usually the euclidean distance.

Several choices for the expression of F can be considered depending on the features to preserve. The simplest one is F (f 0 , .)=f 0 . In the context of image processing, an important feature vector F is provided by images patches i.e.

F (f 0 , u)=F τ (f 0 , u)={f 0 (v) : v∈N τ (u) ∪ {u}}.
In the case of a grayscale image F τ (f 0 , .) is a vector of size (2τ +1) 2 corresponding to the values of f 0 in a square window of size (2τ +1)×(2τ +1) centered at vertex u (a pixel). Color images can be handled using features of dimension 3×(2τ +1) 2 . Then, the resultant weight function directly incorporates local or nonlocal features [START_REF] Elmoataz | Unifying local and nonlocal processing with partial difference operators on weighted graphs[END_REF]. This feature vector has been proposed in the context of texture synthesis [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF], and further used in the context of image processing [START_REF] Buades | Nonlocal image and movie denoising[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF].

Graph Based Discrete Gradients

Let G=(V, E, w) be a weighted graph. The discrete weighted gradient of a function f ∈H(V ) at a vertex u∈V is defined by

(∇ w f )(u) = (∂ v f (u)) (u,v)∈E where ∂ v f (u)= w(u, v)(f (v)-f (u)
) corresponds to the discrete (partial) derivative of f with respect to the edge (u, v). These definitions have been used by [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF] for image and mesh regularization. Based on the latter works, two discrete formulations of weighted morphological gradients on graphs have been proposed by [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF]: namely, the weighted external ∇ + w and the internal ∇ - w gradient operators. For u∈V

(∇ + w f )(u) = (∂ + v f (u)) (u,v)∈E and (∇ - w f )(u) = (∂ - v f (u)) (u,v)∈E , (4)
where the external

∂ + v f (u) and the internal ∂ - v f (u) discrete partial derivatives are ∂ + v f (u) = max(0, ∂ v f (u)) and ∂ - v f (u) = -min(0, ∂ v f (u)) , with ∂ - v f (u)=∂ + u f (v).
When the weight is constant (w=g 0 ) these definitions recover the classical directional derivative operators.

The L p -norm (with 0<p<+∞) and the L ∞ -norm of gradients ( 4) are

(∇ ± w f )(u) p = v∼u w(u, v) p/2 |(f (v)-f (u)) ± | p 1/p and (5)
(∇ ± w f )(u) ∞ = max v∼u w(u, v) 1/2 |(f (v)-f (u)) ± | . (6)
Notation v∼u means that vertex v is adjacent to u. ∇ ± w refers to both external and internal gradient (with respect to the sign) and (a) + = max(0, a) and (a) -= min(0, a). These gradients have the following property:

(∇ w f )(u) p p = (∇ + w f )(u) p p + (∇ - w f )(u) p p with 0<p< + ∞.
Moreover, with a constant weight function and p=∞, (6) recovers the usual expression of algebraic morphological external and internal gradients.

Associated Dilation and Erosion Processes. Continuous morphology (see [START_REF] Brockett | Evolution equations for continuous-scale morphological filtering[END_REF] and references therein) defines flat dilation δ : IR n → IR n and erosion ε : IR n → IR n of a given function f :IR n →IR by a structuring set B={x∈IR n : x p ≤1} with the following general partial differential equations (PDEs):

∂δ(f ) ∂t = ∂ t f = + ∇f p and ∂ε(f ) ∂t = ∂ t f = -∇f p . (7)
Ta et al. [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF] has proposed a discrete version over graphs of these equations by using morphological gradients (4). Given a graph G=(V, E, w), a function f ∈H(V ) and for all u∈V , discrete analogue of [START_REF] Abgrall | Numerical discretization of the first-order hamilton-jacobi equations on triangular meshes[END_REF] are

∂δ(f (u)) ∂t =∂ t f (u)=+ (∇ + w f )(u) p and ∂ε(f (u)) ∂t =∂ t f (u)=-(∇ - w f )(u) p . (8)
Equations ( 8) correspond to dilation and erosion over weighted graphs. They constitute a morphological framework [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF] based on PdEs that extends algebraic and continuous morphological operators for images and high dimensional data processing.

Eikonal Equation on Weighted Graphs

In this Section, we present our formulation to approximate the eikonal equation (3) over weighted graphs by considering PdEs and the morphological gradients presented in the previous Section.

With morphological processes described by [START_REF] Shu | High order WENO schemes for hamiltonjacobi equations on triangular meshes[END_REF], the time dependent eikonal formulation (3) can be viewed as an erosion process regarding the minus sign and a null potential function P . With the corresponding internal gradient (∇ - w ) involved in discrete PdEs based erosion process, (3) can be directly rewritten with weighted graphs. Given a graph G=(V, E, w) and a function f ∈H(V ), we obtain a discrete PdEs based version of the system (3)

     ∂f (u, t)/∂t = -∇ - w f (u) p + P (u) u∈V f (u, t)=φ(u) u∈V 0 ⊂ V f (u, 0)=φ 0 (u) u∈V
, where V 0 corresponds to the initial seed vertices. With f n (u) ≈ f (u, n∆t), this iterative numerical scheme is obtained for all u∈V :

f n+1 (u)=f n (u) -∆t (∇ - w f n )(u) p -P (u) . (9)
The steady state (i.e. given a fixed number n of iteration or when f n+1 -f n < ǫ) of this process is the solution of the eikonal equation [START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF].

Injecting the corresponding internal gradient norm in [START_REF] Mémoli | Fast computation of weighted distance functions and geodesics on implicit hyper-surfaces[END_REF], we obtain for the L p -norm (5) and the L ∞ -norm (6)

f n+1 (u)=f n (u) -∆t v∼u w(u, v) p/2 | min(0, f (v)-f (u))| p 1/p -P (u) , (10)
f n+1 (u)=f n (u) -∆t max v∼u w(u, v) 1/2 | min(0, f (v)-f (u))| -P (u) . (11)
The proposed methodology leads to a simple and common formulation that constitutes an adaptative framework for the eikonal equation. Indeed, our approach only depends on the p value and the weight function w. In Sect. [START_REF] Maragos | Curve evolution, differential morphology and distance transforms as applied to multiscale and eikonal problems[END_REF], experiments show how the framework can be adapted to address image segmentation or data clustering.

Relations with other schemes. Scheme (9) has the advantage to work on any graph structures. Then, with an adapted graph topology and an appropriated weight function, the proposed formulation is linked to well-known schemes such as Osher-Sethian Hamiltonian discretization scheme or the graph based Dikjstra algorithm.

Osher-Sethian scheme. Let G 0 =(V, E, g 0 ) be an unweighted 4-adjacency grid graph associated with an image. Then, [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] recovers the exact Osher-Sethian upwind first order Hamiltonian discretization scheme [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] when p=2 and using G 0 :

f n+1 (u)=f n (u) -∆t v∼u | min(0, f (v)-f (u))| 2 1/2 -P (u) .
Replacing vertices u∈V and their neighborhood by their spatial coordinates (x, y), the latter expression can be rewritten as

f n+1 ((x, y))=f n ((x, y))-∆t | min 0, f n ((x, y))-f n ((x-1, y)) 2 +| max 0, f n ((x+1, y))-f n ((x, y)) 2 +| min 0, f n ((x, y))-f n ((x, y-1)) 2 +| max 0, f n ((x, y+1))-f n ((x, y)) 2 1/2 -P ((x, y)) ,
since min(0, a-b) 2 = max(0, b-a) 2 . This equation corresponds to the discretization scheme of the Hamilton-Jacobi equations proposed by [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. Dikjstra scheme. Let G=(V, E, g 0 ) be an unweighted graph. Then, [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] corresponds to an iterative version of the Dikjstra shortest path algorithm defined on graphs of arbitrary structure. Indeed, in the case where p=∞, ∆t=1 and with G, [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] becomes, for all u∈V

f n+1 (u) = f n (u) -max v∼u | min(0, f (v)-f (u))| + P (u) = min v∼u (f n (v)) + P (u) ,
by considering the neighborhood of u as N (u)∪{u} and with the properties that max(0, a-b)= -min(0, b-a) and min(0, a-b)= min(a, b)-b. This equation corresponds to a shortest path algorithm for a given graph where at each step, the distance f (u) at vertex u corresponds to the minimal distance in its neighborhood.

Experiments

The proposed formulation of the eikonal equation and can be used to process any function defined on vertices of a graph or on any arbitrary discrete domain. This Section illustrates the potentialities of our formulation through examples of weighted distance computation, image segmentation and unorganized high dimensional data processing. Different graph structures and weight functions are also used to show the flexibility of our approach. In the sequel, all experiments are obtained with a constant potential function P =1. Clearly, a different potential function can be adapted for a particular application. The objective of the following experiments is not to solve a particular application. They only illustrate the potential and the behavior of our eikonal equation formulation.

Adaptative Front Propagation and Weighted Distances. Figure 1 shows the adaptivity of our formulation in order to compute weighted distances. Indeed, this example shows results for different p values, graph topologies, weight functions and features F . The initial seed is located at the top left corner of the original grayscale image f 0 :H(V )→IR. First, second and third rows of Fig. 1 show results for p=2, 1 and ∞ respectively, where [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] and ( 11) are used. All the results correspond to color distance maps (red for small and blue for large distances) where iso-levels sets are superimposed in white. First and second columns of Fig. 1 show results obtained with unweighted (w=g 0 ) graphs. First column uses a 4-adjacency grid graph (G 0 ) and corresponds to the classical case. Second column uses a 25-adjacency grid graph (G 2 ) and shows the effect of a larger neighborhood. Third and fourth columns show results obtained with weighted graphs. Third column considers graph G 0 weighted by function g 2 with F =f 0 . By using non constant weights, image information is automatically integrated in the distance computation that modifies the front evolution speed particularly into the textured sub-image. Fourth column shows the nonlocal case where graph G 2 is constructed and weighted with function g 2 associated with patches of size 11×11. In that case, repetitive information are clearly captured by the weights that stops the front propagation around the textured sub-image. Finally, segmentation of the textured sub-image can be simply obtained by thresholding the computed distances.

Image Segmentation with Region Based Graphs. The goal of the following two examples is not to show a perfect segmentation but to show how we can take advantage of graph topologies in image segmentation. The basic idea is to consider that image pixels are not the only relevant components in image Figures 2(c) and 3(c) are reconstructed images from partitions with the mean color value for each region. Figure 2 presents an example of image segmentation based on RAG and also shows that this graph structure can accelerate segmentation processes. This example compares segmentation obtained by a 4-adjacency grid graph G 0 weighted by function g 2 with pixel grayscale values (Fig. 2(d)) and segmentation result with a RAG constructed from partition (Fig. 2(b)) and weighted by function g 2 with mean values (Fig. 2(e)). Color distance maps are obtained with the initial seeds (white points) in Fig. 2(a). Segmentations are performed by thresholding the obtained distances. Results show similar behaviors both on distance maps and segmentations while drastically speeding-up the segmentation process in the RAG case. Indeed, the number of vertices in the RAG represents approximatively 3% as compared to the number of vertices in the pixel based graph. The direct consequence is a decreasing of the computational complexity thanks to the reduced amount of data to consider. On a standard computer the computing time can be decreased by a 10 factor.

Figure 3 shows another benefit of using a RAG structure: nonlocal (non spatially connected) object segmentation. This experiment compares segmentation results with RAG (Fig. 3(d)) and nonlocal RAG (Fig. 3 In nonlocal RAG case, each vertex neighborhood is extended by a 5 nearest neighborhood based on mean value feature. The obtained graph is a RAG∪5-NNG graph. Figures 3(d) and 3(e) show color distance maps computed from initial seeds (white stroke) in Fig. 3(a) and final segmentations. For local case (Fig. 3(d)), object marked by seeds is well segmented with respect to close distances (red color). The other objects are far (blue color) and the final segmentation only extracts the marked one. For nonlocal case (Fig. 3(e)), the distance within the marked object is close to the initial seeds. In addition the distances to other triangles in the scene are also computed as close to seeds (red color). The consequence is that all the objects in the image are extracted by thresholding even if they are not spatially close with a minimal number of initial seeds.

Unorganized High Dimensional Data Processing The following experiments show applications of our formulation of the eikonal equation for the processing of high dimensional data in irregular domains. Figure 4 shows applications of the eikonal equation for data clustering and shortest path problems. The initial data set (Fig. 4(a)) is constituted of 133 images of head pose. Each image is of size 29×29. From this data set, two possible applications can be performed: clustering and head pose transition estimation. The goal here is not to solve machine leaning problems, but to show that these problems can be addressed by our formulation of eikonal equation. In order to process such data, a graph (|V |=133) is constructed where each vertex represents an image and is described by a feature of size 29×29 (i.e IR 841 ). In the following results, initial seeds (images) are represented with white boundaries. Points that are close and far to seeds are respectively represented with blue and red colors in distance maps (Fig. 4(b) and 4(c)). Figures 4(b) and 4(d) show the application of the eikonal equation for data clustering. Such an application can be used for data set exploration or semisupervised learning: given an input seed (query) one wants to obtained the closest points with respect to the initial input. Figures 4(c) and 4(e) shows another example of application of the eikonal equation for data set. Given two initial images, one wants to recover a transition sequence of images that separates them. This problem can be viewed as a shortest path problem solved by the eikonal equation. Figure 4(c) shows the distance map obtained from the initial seeds. Figure 4(e) shows the obtained path from seed at top left to seed at bottom right.

These experiments show satisfying results and the ability of our approach to address machine learning problems even if a simple euclidean distance is used to compare data points. Clearly, results can be improved by using well adapted distances or features estimation.

Conclusion

In this paper, a discrete version of the eikonal equation over weighted graphs of arbitrary structure is proposed. Solution of the eikonal equation based on PdEs, discrete gradients and weighted graphs is presented. The proposed formulation constitutes a simple, common and adaptative framework that recovers well-known definitions and unifies local and nonlocal configurations in the context of image processing. This framework can consider any discrete data that can be represented by weighted graphs. Through experiments, we have shown the potentiality and the flexibility of our approach to address image segmentation and unorganized high dimensional data processing. Finally, an ongoing work is to address the stationary (time independent) version of the eikonal equation and to solve this equation by considering fast marching like methods on arbitrary graphs within our framework.

Fig. 1 .

 1 Fig. 1. Front propagation and weighted distances with different p values, graph configurations G, weights w and features F . Figures represent color distance maps with iso-level sets obtained by thresholding the distances. The seed is located at the top left corner (see text for more details).

  (e)). Both graphs are computed from partition 3(b) weighted by function g 2 with mean color values.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Acceleration of image segmentation process. (a) original image (150×235) with 35 250 pixels. (b) partition with 999 regions (97% of reduction in terms of image components). (c) reconstructed image with mean color value. (d) and (e): at left, distance color maps (red for small and blue for large distances) and at right, final segmentations. Images (d) are obtained with a pixel based graph computed from (a). Images (e) are obtained with a RAG constructed with (b) and (c) (see text for more details).

Figure 4 (

 4 b) shows the distance map obtained from a single initial seed. Figure 4(d) shows clustering results. Initial input has a white boundary. The 10 closest images are located at the top and the 10 farthest are located at the bottom of Fig. 4(d).

Fig. 4 .

 4 Fig. 4. High dimensional data clustering and shortest path. (b) and (c) color distance maps (blue for small and red for large distance) images superimposed with white boundaries are initial seeds. (d) clustering results where at top the 10 closest and at bottom the 10 farthest with respect to the seed (white boundary). (e) shortest path from the two initial seeds (white boundary).