Pawel Mlynarski

Hervé Delingette

Antonio Criminisi

Nicholas Ayache

Deep Learning with Mixed Supervision for Brain Tumor Segmentation

Keywords: Semi-supervised learning, Convolutional Neural Networks, segmentation, tumor, MRI

Most of the current state-of-the-art methods for tumor segmentation are based on machine learning models trained on manually segmented images. This type of training data is particularly costly, as manual delineation of tumors is not only time-consuming but also requires medical expertise. On the other hand, images with a provided global label (indicating presence or absence of a tumor) are less informative but can be obtained at a substantially lower cost. In this paper, we propose to use both types of training data (fully-annotated and weakly-annotated) to train a deep learning model for segmentation. The idea of our approach is to extend segmentation networks with an additional branch performing image-level classification. The model is jointly trained for segmentation and classification tasks in order to exploit information contained in weakly-annotated images while preventing the network to learn features which are irrelevant for the segmentation task. We evaluate our method on the challenging task of brain tumor segmentation in Magnetic Resonance images from BRATS 2018 challenge. We show that the proposed approach provides a significant improvement of segmentation performance compared to the standard supervised learning. The observed improvement is proportional to the ratio between weakly-annotated and fully-annotated images available for training.

Introduction

Cancer is today the third cause of mortality worldwide. In this paper, we focus on segmentation of gliomas, which are the most frequent primary brain cancers [START_REF] Goodenberger | Genetics of adult glioma[END_REF]. Gliomas are particularly malignant tumors and can be broadly classified according to their grade into low grade gliomas (grades I and II defined by World Health Organization) and high grades gliomas (grades III-IV). Glioblastoma multiforme is the most malignant form of glioma and is associated with a very poor prognosis: the average survival time under therapy is between 12 and 14 months.

Medical images play a key role in diagnosis, therapy planning and monitoring of cancers. Treatment protocols often include evaluation of tumor volumes and locations. In particular, for radiotherapy planning, clinicians have to manually delineate target volumes, which is a difficult and timeconsuming task. Magnetic Resonance (MR) images [START_REF] Bauer | A survey of mri-based medical image analysis for brain tumor studies[END_REF] are particularly suitable for brain cancer imaging. Different MR sequences (T2, T2-FLAIR, T1, T1+gadolinium) highlight different tumor subcomponents such as edema, necrosis or contrast-enhancing core.

In recent years, machine learning methods have achieved impressive performance in a large variety of image recognition tasks. Most of the recent state-of-the-art segmentation methods are based on Convolutional Neural Networks (CNN) [START_REF] Lecun | Convolutional networks for images, speech, and time series, The handbook of brain theory and neural networks[END_REF][START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. CNNs have the considerable advantage of automatically learning relevant image features. This ability is particularly important for the tumor segmentation task. CNN-based methods [START_REF] Pereira | Deep convolutional neural networks for the segmentation of gliomas in multi-sequence mri[END_REF][START_REF] Kamnitsas | Efficient multi-scale 3d cnn with fully connected crf for accurate brain lesion segmentation[END_REF][START_REF] Kamnitsas | Ensembles of multiple models and architectures for robust brain tumour segmentation[END_REF][START_REF] Wang | Automatic brain tumor segmentation using cascaded anisotropic convolutional neural networks[END_REF] have obtained the best performances on the four last editions of Multimodal Brain Tumor Segmentation Challenge (BRATS) [START_REF] Menze | The multimodal brain tumor image segmentation benchmark (brats)[END_REF][START_REF] Bakas | Advancing the cancer genome atlas glioma mri collections with expert segmentation labels and radiomic features[END_REF].

Most of the segmentation methods based on machine learning rely uniquely on manually segmented images. The cost of this annotation is particularly high in medical imaging where manual segmentation is not only timeconsuming but also requires high medical competences. Image intensity of cancerous tissues in MRI or CT scans is often similar to the one of surrounding healthy or pathological tissues, making the exact tumor delineation difficult and subjective. In the case of brain tumors, according to [START_REF] Menze | The multimodal brain tumor image segmentation benchmark (brats)[END_REF], the interrater overlap of expert segmentations is between 0.74 and 0.85 in terms of Dice coefficient. For these reasons, high-quality manual tumor segmentations are generally available in very limited numbers. Segmentation approaches able to exploit images with weaker forms of annotations are therefore of particular interest.

In this paper, we assume that the training dataset contains two types of images: fully-annotated (with provided ground truth segmentation) and weakly-annotated, with an image-level label indicating presence or absence of a tumor tissue within the image (Fig. 1). We refer to this setting as 'mixed supervision'. The latter type of annotations can be obtained at a substantially lower cost as it is less time-consuming, potentially requires less medical expertise and can be obtained without the use of a dedicated software.

We introduce a novel CNN-based segmentation model which can be trained using weakly-annotated images in addition to fully-annotated images. We propose to extend segmentation networks, such as U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], with an additional branch, performing image-level classification. The model is trained jointly for both tasks, on fully-annotated and weakly-annotated images. The goal is to exploit the representation learning ability of CNNs to learn from weakly-annotated images while supervising the training using fully-annotated images in order to learn features relevant for the segmentation task. Our approach differs from the standard semi-supervised learning as we consider weakly-annotated data instead of totally unlabelled data. To the best of our knowledge, we are the first to combine pixel-level and image-level labels for training of models for tumor segmentation.

We perform a series of cross-validated tests on the challenging task of segmentation of gliomas in MR images from BRATS 2018 challenge. We evaluate our model both for binary and multiclass segmentation using a variable number of ground truth segmentations available for training. Since all 3D images from the BRATS 2018 contain brain tumors, we focus on the 2D problem of tumor segmentation in axial slices of a MRI and we assume slice-level labels for weakly-annotated images. Using approximately 220 MRI with slice-level labels and a varying number of fully-annotated MRI, we show that our approach significantly improves the segmentation accuracy when the number of fully-annotated cases is limited.

Related work

In the literature, there are several works related to weakly-supervised and semi-supervised learning for object segmentation or detection. Most of the related works were applied to natural images.

The first group of weakly-supervised methods aims to localize objects using only weakly-annotated images for training. When only image-level labels are available, one approach is to design a neural network which outputs two feature maps per class (interpreted as 'heat maps' of the class) which are then pooled to obtain an image-level classification score penalized during the training [START_REF] Pathak | Fully convolutional multiclass multiple instance learning[END_REF][START_REF] Pinheiro | From image-level to pixel-level labeling with convolutional networks[END_REF][START_REF] Saleh | Built-in foreground/background prior for weakly-supervised semantic segmentation[END_REF][START_REF] Bearman | Whats the point: Semantic segmentation with point supervision[END_REF][START_REF] Wang | Automated detection of clinically significant prostate cancer in mp-mri images based on an end-to-end deep neural network[END_REF]. At test time, these 'heat maps' are used for detection (determining a bounding box of the object) or segmentation. To guide the training process, some works use self-generated spatial priors [START_REF] Pinheiro | From image-level to pixel-level labeling with convolutional networks[END_REF][START_REF] Saleh | Built-in foreground/background prior for weakly-supervised semantic segmentation[END_REF][START_REF] Bearman | Whats the point: Semantic segmentation with point supervision[END_REF] or inconsistency measures [START_REF] Wang | Automated detection of clinically significant prostate cancer in mp-mri images based on an end-to-end deep neural network[END_REF] in the loss function. To obtain an image-level score, in [START_REF] Pathak | Fully convolutional multiclass multiple instance learning[END_REF][START_REF] Bearman | Whats the point: Semantic segmentation with point supervision[END_REF], global maximum pooling is used. Application of the maximum function on large feature maps may cause optimization problems as training of neural networks is based on computation of gradients [START_REF] Dreyfus | Artificial neural networks, back propagation, and the kelley-bryson gradient procedure[END_REF]. LogSumExp approximation of the maximum [START_REF] Boyd | Convex optimization[END_REF] is therefore used in the works [START_REF] Pinheiro | From image-level to pixel-level labeling with convolutional networks[END_REF][START_REF] Saleh | Built-in foreground/background prior for weakly-supervised semantic segmentation[END_REF] in order to partially limit this problem. Average pooling on small feature maps was used by Wang et al [START_REF] Wang | Automated detection of clinically significant prostate cancer in mp-mri images based on an end-to-end deep neural network[END_REF] for the problem of detection of lung nodules.

Another type of weakly-supervised methods aims to detect objects in natural images based on classification of image subregions [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF][START_REF] Oquab | Is object localization for free?-weakly-supervised learning with convolutional neural networks[END_REF] using pretrained classification networks such as VGG-Net [START_REF] Simonyan | Very deep convolutional networks for largescale image recognition[END_REF] or AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. In fact, one particularity of natural images is their recursive aspect: one image can correspond to a subpart of another image (e.g. two images of the same object taken from different distances). A classification network trained on a large dataset may therefore be used on a subregion of a new image in order to determine if it contains an object of interest.

Pre-trained classification networks were also used to detect objects by determining image subregions whose modification influences the global clas-sification score of a class. In [START_REF] Simonyan | Deep inside convolutional networks: Visualising image classification models and saliency maps[END_REF], Simonyan et al. propose to compute the gradient of the classification score with respect to the intensities of pixels and to threshold it in order to localize the object of interest. However, these partial derivatives represent a very weak information for tumor segmentation, which requires a complex analysis of the spatial context. The method proposed in [START_REF] Bergamo | Self-taught object localization with deep networks[END_REF] is based on replacing image subregions by the mean value in order to measure the drop of the classification score.

Overall, the reported segmentation performances of weakly-supervised methods are considerably lower than the ones obtained by semi-supervised and supervised approaches. In absence of pixel-level labels, a model may learn irrelevant features, due for example to co-occurrences of objects or image acquisition differences in the case of multicenter medical data. Despite the cost of manual segmentation, at least few fully-annotated images can still be obtained in many cases.

In standard semi-supervised learning [START_REF] Cheplygina | Not-so-supervised: a survey of semi-supervised, multi-instance, and transfer learning in medical image analysis[END_REF] for classification, the training data is composed both of labelled samples and unlabelled samples. Unlabelled samples can be used to encourage the model to satisfy some properties on relations between labels and the feature space. Common properties include smoothness (points close in the feature space should be close in the target space), clustering (labels form clusters in the feature space) and low density separation (decision boundaries should be in low density regions of the feature space). Semi-supervised learning based on these properties can be performed by graph-based methods such as the recent work of Kamnitsas et al. [START_REF] Kamnitsas | Semi-supervised learning via compact latent space clustering[END_REF]. The main idea of such methods is to propagate labels in a fully-connected graph whose nodes are samples (labelled and unlabelled) and whose edges are weighted by similarities between samples. The use of graph-based semi-supervised methods is difficult for segmentation, in particular because it implies computation of similarity metrics between samples, whereas each single image is generally composed of millions of samples (pixels or voxels).

Relatively few works were proposed for semi-supervised learning for image segmentation. Some semi-supervised approaches are based on self-training, i.e. training of a machine learning model on self-generated labels. Iterative algorithms similar to EM [START_REF] Zhang | Segmentation of brain mr images through a hidden markov random field model and the expectation-maximization algorithm[END_REF] were proposed for natural images [START_REF] Papandreou | Weakly-and semi-supervised learning of a deep convolutional network for semantic image segmentation[END_REF] and medical images [START_REF] Rajchl | Deepcut: Object segmentation from bounding box annotations using convolutional neural networks[END_REF]. Recently, Hung et al [START_REF] Hung | Adversarial learning for semi-supervised semantic segmentation[END_REF] proposed a method based on Generative Adversarial Networks [START_REF] Goodfellow | Generative adversarial nets[END_REF] where the generator network performs image segmentation and the discriminator network tries to determine if a segmentation corresponds to the ground truth or the segmentation produced by the generator. The discriminator network is used to produce confidence maps for self-training. The approaches based on self-training have the drawback of learning on uncertain labels (produced by the model itself) and training of such models is difficult.

Other approaches assume mixed levels of supervision similarly to our approach. Hong et al [START_REF] Hong | Decoupled deep neural network for semisupervised semantic segmentation[END_REF][START_REF] Hong | Learning transferrable knowledge for semantic segmentation with deep convolutional neural network[END_REF] proposed decoupled classification and segmentation, an approach for segmentation of objects in natural images based on a two-step training with a varying level of supervision. This architecture is composed of two separate networks trained sequentially, one performing image-level classification and used as encoder, and the another one taking as input small feature maps extracted from the encoder and performing segmentation. An important drawback of such design, in the case of tumor segmentation, is that the segmentation network does not take as input the original image and can therefore miss important details of the image (e.g. small tumors).

Our approach is related to multi-task learning [START_REF] Evgeniou | Regularized multi-task learning[END_REF]. In our case, the goal of training for two tasks (segmentation and classification) is to exploit all the available labels and to guide the training process to learn relevant features. The approach closest to ours is the one of Shah et al. [START_REF] Shah | Ms-net: Mixed-supervision fullyconvolutional networks for full-resolution segmentation[END_REF]. In this work, the authors consider three types of annotations: segmentations, bounding boxes and seed points at the borders of objects. A neural network is trained using these three types of training data. In our work, we exploit the use of a significantly weaker form of annotations, image-level labels.

Joint classification and segmentation with Convolutional Neural Networks

Deep learning model for binary segmentation

We designed a novel deep learning model, which aims to take advantage of all available voxelwise and image-level annotations. We propose to extend a segmentation CNN with an additional subnetwork performing image-level classification and to train the model for the two tasks jointly. Most of the layers are shared between the classification and segmentation subnetworks in order to transfer the information between the two subnetworks. In this paper we present the 2D version of our model, which can be used on different types of medical images such as slices of a CT scan or a multisequence MRI.

The proposed network takes as input an image of dimensions 300x300 and extends U-net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] which is currently one of the most used architectures for segmentation tasks in medical imaging. This segmentation architecture is composed of an encoder part and a decoder part which are connected by concatenations between layers at the same scale, in order to combine lowlevel and local features with high-level and global features. This design is well suited for the tumor segmentation task since the classification of a voxel as tumor requires to compare its value with its close neighborhood but also taking into account a large spatial context. The last convolutional layer of Unet produces pixelwise classification scores, which are normalized by softmax function during the training phase. We apply batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] in all convolutional layers except the final layer.

We propose to add an additional branch to the network, performing image-level classification (Fig. 2), in order to exploit the information contained in weakly-annotated images during the training. This classification branch takes as input the second to last convolutional layer of U-net (representing a rich information extracted from a local and a long-range spatial context) and is composed of one mean-pooling, one convolutional layer and 7 fully-connected layers.

The goals of taking a layer from the final part of U-Net as input of the classification branch are both to guide the image-level classification task and to force the major part of the segmentation network to take into account weakly-annotated images. This also helps the optimization process by taking advantage of the connectivity of layers in U-Net, helping the flow of gradients of the loss function during the training (in particular, note the connection between the first part and the last part of U-Net).

The second to last layer of the segmentation network outputs 64 feature maps of size 101x101 from which the classification branch has to output two global (image-level) classification scores (tumor absent/tumor present). We first reduce the size of these feature maps by applying a mean-pooling with kernels of size 8x8 and the stride of 8x8. We use the mean pooling rather than max-pooling in order to avoid information loss and optimization problems. One convolutional layer, with ReLU activation and kernels of size 3x3, is then added to reduce the number of feature maps from 64 to 32. The resulting 32 feature maps of size 11x11 are the input of the first fully-connected layer of the classification branch.

According to our experiments, a relatively deep architecture of the classification branch with a limited number of parameters and a skip-connection between layers yields the best performance. This observation is in agreement with current common designs of neural networks. Deep networks have the capacity to learn more complex features, due to applied non-linearities. The connectivity between layers at different depths helps the optimization process (e.g. Res-Net [START_REF] He | Deep residual learning for image recognition[END_REF]). In our case, we use 7 fully-connected layers with ReLU activations (except the final layer) and we concatenate the outputs of the first and the fifth fully-connected layer. The last fully-connected layer outputs image-level classification scores (tumor tissue absent or present).

The model is trained both on fully-annotated and weakly-annotated images for the two tasks jointly (segmentation and classification). We can distinguish between three types of training images. First, images containing a tumor and with provided ground truth segmentation are the most costly ones. The second type are images that do not contain tumor, which implies that none of their pixels corresponds to a tumor. In this case, the ground truth segmentation is simply the zero matrix. The only problematic case is the third one, when the image is labelled as containing a tumor but without provided segmentation.

To train our model, we propose to form training batches containing the three mentioned types of images: k positive cases (containing a tumor) with provided segmentation, m negative cases and n positive cases without provided segmentation. In our experiments we chose k = 4, m = 2, n = 4.

Given a training batch b and the network parameters θ, we use a weighted pixelwise cross-entropy loss on images of types 1 and 2:

Loss b s (θ) = - 1 P k+m i=1 (x,y) w i (x,y) log(p l i,(x,y) (θ)) ( 1
)
where P is the number of pixels; p l i,(x,y) is the classification score given by the network to the ground truth label for pixel (x,y) of the i th image of the batch and w i (x,y) is the weight given to this pixel. The weights are used to limit the effect of class imbalance, since tumor pixels represent a small portion of the image. Weights of pixels are set automatically according to the composition of the training batch (number of pixels of each class) so that pixels associated with healthy tissues have a total weight of t 0 in the loss function and the pixels of the tumor class have a total weight of t 1 , where t 0 and t 1 are target weights fixed manually (t 0 = 0.7 and t 1 = 0.3 in our experiments). It means that if the training batch contains N t pixels labelled as tumor, then each tumor pixel has a weight of t 1 /N 1 (the pixelwise weight is high when the number of tumor pixels is low). We fix a higher target weight for healthy tissues to avoid oversegmentation, while increasing the relative weight of tumor pixels (t 1 = 0.3) compared to the standard non-weighted cross-entropy loss where the relative weight of the tumor class is proportional to the number of tumor pixels, i.e approximately 1-2 %. This type of loss function was used in our previous work [START_REF] Mlynarski | 3d convolutional neural networks for tumor segmentation using long-range 2d context[END_REF].

The classification loss is a standard cross-entropy loss on all images of the training batch: Loss c = -1 k+m+n k+m+n i=1 log(p l i (θ)) where p l i is the global classification score given by the network to the ground truth global label for the i th image of the batch. In particular, fully-annotated images are also used for training of the classification branch in order to transfer the knowledge from the segmentation task to the image-level classification. We do not apply weights on the classification loss as image-level labels are balanced through the sampling of training batches (having a fixed number of nontumor images).

Since both segmentation and classification losses are normalized, we define the total loss as a convex combination of the classification and segmentation losses: Loss = a * Loss s + (1 -a) * Loss c . In our experiments we fixed a=0.7, i.e. we give a higher importance to segmentation errors. The determination of this parameter was performed empirically through a series of cross-validated tests.

We train our model with a variant of Stochastic Gradient Descent (SGD) with momentum [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF], used also in our previous work [START_REF] Mlynarski | 3d convolutional neural networks for tumor segmentation using long-range 2d context[END_REF]. The main differences with the standard SGD are to divide the gradient by its norm and to compute gradients on several training batches in each iteration, in order to take into account many training examples while bypassing GPU memory constraints.

Extension to the multiclass problem

We extend our model to the multiclass case where each pixel has to be labelled with one of K classes, such as the four ones considered in BRATS challenge (non-tumor, contrast-enhancing core, edema, non-enhancing core). We now assume that image-level labels are provided for each class (absent/present in the image).

Extension of the segmentation subnetwork to the multiclass problem is straightforward, by changing the number of final feature maps to match the number of classes. However, image-level labels are not exclusive, i.e. an image may contain several tumor subclasses. For this reason, we propose to consider one image-level classification output per tumor subclass, indicating absence or presence of the given subclass.

According to our experiments, better performances are obtained when each subclass has its dedicated entire classification branch (Fig. 3). A possible reason is that the image-level classification of tumor subclasses is a challenging task requiring a sufficient number of dedicated parameters.

Training batches are sampled similarly to the binary case, however each tumor subclass has to be present at least once in each training batch. In our implementation, we store lists of paths of images containing tumor subclasses in order to sample from these lists during the training of the model.

In the segmentation loss we empirically fix the following target weights for the four classes (non-tumor, non-enhancing tumor core, edema, enhancingcore): t 0 = 0.7, t 1 = 0.1, t 2 = 0.1, t 3 = 0.1 (all tumor subclasses have an equal weight in the loss function). The loss associated with each classification branch is the same as in the binary case and the total classification loss is the average across all classification branches. We observe the need to decrease the parameter a (weight of the segmentation loss) compared to the binary case. We fixed a = 0.3 according to performed cross-validated experiments.

Experiments

Data

We evaluate our method on the challenging task of brain tumor segmentation in multisequence MR scans, using the Training dataset of BRATS 2018 challenge. It contains 285 multisequence MRI of patients diagnosed with lowgrade gliomas or high-grade gliomas. For each patient, manual ground truth segmentation is provided. In each case, four MR sequences are available: T1, T1+gadolinium, T2 and FLAIR (Fluid Attenuated Inversion Recovery). Preprocessing performed by the organizers includes skull-stripping, resampling to 1 mm 3 resolution and registration of images to a common brain atlas. The resulting volumes are of size 240x240x155. The images were acquired in 19 different imaging centers. In order to normalize image intensites, each image is divided by the median of non-zero voxels (which is supposed to be less affected by the tumor zone than the mean) and multiplied the image by a fixed constant.

Each voxel is labelled with one of the following classes: non-tumor (class 0), contrast-enhancing core (class 3), non-enhancing core (class 1), edema (class 2). The benchmark of the challenge groups classes in three regions: Given that all 3D images of the database contain tumors (no negative cases to train a 3D classification network), we consider the 2D problem of tumor segmentation in axial slices of the brain.

Test setting

The goal of our experiments is to compare our approach with the standard supervised learning. In each of the performed tests, our model is trained on fully-annotated and weakly-annotated images and is compared with the standard U-Net trained on fully-annotated images only. The goal is to compare our model with a commonly used segmentation model on a publicly available database.

We consider three different training scenarios, with a varying number of patients for which we assume a provided manual tumor segmentation. In each scenario we perform a 5-fold cross-validation. In each fold, 57 patients are used for test and 228 patients are used for training. Among the 228 training images, few cases are assumed to be fully-annotated and the remaining ones are considered to be weakly-annotated, with slice-level labels. The fullyannotated images are different in each fold. If the 3D volumes are numbered from 0 to 284, then in k th fold, the test images correspond to the interval [(k-1)*57, k*57 -1], the next few images correspond to fully-annotated images and the remaining ones are considered as weakly-annotated (the folds are generated in a circular way). In the following, FA denotes the number of fully-annotated cases and WA denotes the number of weakly-annotated cases (with slice-level labels).

In the first training scenario, 5 patients are assumed to be provided with a manual segmentation and 223 patients have slice-level labels. In the second and the third scenario, the numbers of fully-annotated cases are respectively 15 and 30 and the numbers of weakly-annotated images are therefore respectively 213 and 198. The three training scenarios are independent, i.e. folds are re-generated randomly (the list of all images is permuted randomly and the folds are generated). In fact, results are likely to depend not only on the number of fully-annotated images but also on qualitative factors (for example the few fully-annotated images may correspond to atypical cases), and the goal is to test the method in various settings. Overall, our approach is compared to the standard supervised learning on 60 tests (5-fold crossvalidation, three independent training scenarios, three binary problems and one multiclass problem).

We evaluate our method both on binary segmentations problems (separately for each of three tumor regions considered in the challenge) and on the end-to-end multiclass segmentation problem. In each binary case, the model is trained for segmentation and classification of one tumor region (whole tumor, tumor core or enhancing core).

Segmentation performance is expressed in terms of Dice score quantifying the overlap between the ground truth (Y ) and the output of a model ( Ỹ ):

DSC( Ỹ , Y ) = 2| Ỹ ∩ Y | | Ỹ | + |Y | (2)

Results

The main observation is that our model with mixed supervision provides a significant improvement over the standard supervised approach (U-Net trained on fully-annotated images) when the number of fully-annotated images is limited. In the two first training scenarios (5 FA and 15 FA), our model outperformed the supervised approach on the three binary segmentation problems (Table 1) and in the multiclass setting (Table 3). The largest 2 for the binary problems and in Table 4 for the multiclass problem. Our approach provided an improvement in all folds of the second scenario and for all tumor regions, except one fold for enhancing core in the binary setting.

In the third scenario (30 FA + 198 WA), our approach and the standard supervised approach obtained similar performances. Segmentation performance increases quickly with the first fully-annotated cases, both for the standard supervised learning and the learning with mixed supervision. For instance, mean Dice score obtained by the supervised approach for whole tumor increases from 70.39, in the case with 5 fully-annotated images, to 77.9 in the case with 15 fully-annotated images. Our approach using 5 fully-annotated images and 223 weakly-annotated images obtained a slightly better performance (78.3) than the supervised approach using 15 fully-annotated cases (77.9). This result is represented on Fig. 7.

Note that each fully-annotated case corresponds to a large 3D volume with voxelwise annotations. Each manually segmented axial slice of size 240x240 corresponds to 57 600 labels, which represents indeed a huge amount In particular, in the first training scenario, the standard supervised approach does not detect the tumor core zone, in contrast to our method.

of information compared to one global label simply indicating presence of absence of a tumor tissue within the slice.

In terms of the annotation cost, manual delineation of tumor tissues in one MRI may take about 45 minutes for an experienced oncologist using a dedicated segmentation tool. Determining the range of axial slices containing tumor tissues may take 1-2 minutes but can be done without a specialized software. More importantly, determining global labels may require less medical expertise than performing an exact tumor delineation and can therefore be performed by a larger community.

Conclusion and future work

In this paper we proposed a new deep learning approach for tumor segmentation which takes advantage of weakly-annotated medical images during the training of neural networks, in addition to a small number of manually segmented images. In our approach, we propose to use neural networks producing both voxelwise and image-level outputs. The classification and segmentation subnetworks share most of their layers and are trained jointly using both fully-annotated and weakly-annotated data. We performed a large number of cross-validated experiments to test our method both in binary and multiclass settings. Our experiments showed that the use of weaklyannotated data improves the segmentation performance significantly when the number of manually segmented images is limited. Our model is end-toend and straightforward to implement with common deep learning libraries such as Theano [START_REF] Bergstra | Theano: A cpu and gpu math compiler in python[END_REF] or TensorFlow [START_REF] Abadi | Tensorflow: Large-scale machine learning on heterogeneous distributed systems[END_REF]. The code of our method will be made publicly available in order to encourage other researchers to continue the research in the field.

An interesting step of the future work would be to extend our method to an end-to-end 3D segmentation. In our paper we focused on the 2D segmentation problem, in particular because all 3D images from the BRATS 2018 database contain tumors whereas we also need non-tumor images to train the classification part of our model. One advantage of the 3D application would be to take into account a richer spatial context in the case of MRI or CT scans. Furthermore, volume-level labels require less effort than slice-level labels and would therefore be easier to obtain, even if these labels are also less informative.

In our tests, we used approximately 220 weakly-annotated MRI, which is relatively a limited number. An important future step would be to test our method on a database containing a considerably larger number of weaklyannotated images (thousands, millions).

Figure 1 :

 1 Figure 1: Different levels of supervision for training of segmentation models. Standard models are trained on fully-annotated images only, with pixel-level labels. Weaklysupervised approaches aim to train models using only weakly-annotated images, e.g. with image-level labels. Our model is trained with a mixed supervision, exploiting both types of training images.

Figure 2 :

 2 Figure 2: Architecture of our model for binary segmentation. The numbers of outputs are specified below boxes representing layers. The height of rectangles represents the scale (increasing with pooling operations). The dashed lines represent concatenation operations. The proposed architecture is an extended version of U-net, with a subnetwork performing image-level classification. Training of the model corresponds to a joint minimization of two loss functions related respectively to segmentation and image-level classifcation tasks.

Figure 3 :

 3 Figure 3: Extension of our model to the multiclass problem. The number of final feature maps of the segmentation subnetwork is equal to the number of classes (4 in our case). As image-level labels (class present/absent) are not exclusive, we consider one classification branch per tumor subclass.

Figure 4 :

 4 Figure 4: Examples of multisequence MRI from the BRATS 2018 database. While T2 and T2-FLAIR highlight the edema induced by the tumor, T1 is suitable for determining the tumor core. In particular, T1 acquired after injection of a contrast product (T1c) highlights the tumor angiogenesis, indicating presence of highly proliferative cancer cells.

Figure 5 :

 5 Figure 5: Comparison of our approach with the standard supervised learning for binary segmentation of the 'whole tumor' region. Each row represents the same test example (first image of Fig. 4) from a different training scenario (5, 15 or 30 fully-annotated scans available for training). FA and WA refer respectively to the number of fully-annotated MRI and weakly-annotated MRI (with slice-level labels). The results are displayed on MRI T2-FLAIR sequence. The performance of both models improves with the number of manual segmentations available for training.

Figure 6 :

 6 Figure 6: Comparison of our approach with the standard supervised learning for binary segmentation of the 'tumor core' region (test example corresponding to the bottom image of Fig. 4). Each row corresponds to a different training scenario (5, 15 or 30 fullyannotated scans available for training). FA and WA refer to the numbers of fully-annotated and weakly-annotated scans. The results are displayed on MRI T1+gadolinium. The observations are similar to the problem of binary segmentation of the 'whole tumor' region.In particular, in the first training scenario, the standard supervised approach does not detect the tumor core zone, in contrast to our method.

Figure 7 :

 7 Figure 7: Illustration of the improvement provided by the mixed supervision for binary segmentation of the 'whole tumor' region. Mixed supervision using 5 fully-annotated MRI and 223 weakly-annotated MRI obtains a slightly better performance than the standard supervised approach using 15 fully-annotated MRI. The improvement provided by the weakly-annotated images decreases with the number of available ground truth segmentations.

Table 1 :

 1 Mean Dice scores (5-fold cross-validation, 57 test cases in each fold) in the three binary segmentation problems obtained by the standard supervised approach and by our model trained with mixed supervision.

		Whole Tumor Tumor Core	Enhancing core
	Standard supervision 5 FA	70.39	48.14	55.74
	Mixed supervision 5 FA + 223 WA	78.34	50.11	60.06
	Standard supervision 15 FA	77.91	58.33	62.88
	Mixed supervision 15 FA + 213 WA	80.92	63.23	66.61
	Standard supervision 30 FA	83.95	66.17	69.15
	Mixed supervision 30 FA + 198 WA	83.84	68.30	67.18
	Standard supervision 228 FA	86.80	77.09	72.20

Table 2 :

 2 Results obtained for the three binary problems (whole tumor, tumor core, enhancing core) on different folds in the case with 15 fully-annotated images and 213 weaklyannotated images. FA) for the whole tumor region where the improvement is of 8 points of the mean Dice score in the binary setting and of 9 points of Dice in the multiclass setting. Results on different folds of the second scenario (intermediate case, 15 FA) are displayed in Table

		fold1	fold2	fold3	fold4	fold5	mean
	Standard supervision, whole tumor	76.23	78.15	78.13	77.67	79.35	77.91
	Mixed supervision, whole tumor	82.36 81.03	78.96 79.88 82.35 80.92
	Standard supervision, tumor core	61.46	61.17	56.68	56.42	55.94	58.33
	Mixed supervision, tumor core	63.15 66.82	63.45 60.83 61.91 63.23
	Standard supervision, enhancing core	66.33	61.08	57.86	68.09	61.02	62.88
	Mixed supervision, enhancing core	68.72 70.65	60.34	67.55	65.80 66.61
	improvements are in the first scenario (5

Table 3 :

 3 Mean Dice scores (5-fold cross-validation, 57 test cases in each fold) obtained by the standard supervised approach and by our model in the multiclass setting.

		Whole Tumor Tumor Core	Enhancing core
	Standard supervision 5 FA	67.61	51.12	58.15
	Mixed supervision 5 FA + 223 WA	76.64	56.30	58.19
	Standard supervision 15 FA	74.46	59.87	61.85
	Mixed supervision 15 FA + 213 WA	79.39	63.91	65.71
	Standard supervision 30 FA	81.10	67.48	68.67
	Mixed supervision 30 FA + 198 WA	81.23	66.33	67.69
	Standard supervision 228 FA	85.67	78.78	74.14

Table 4 :

 4 Results obtained in the multiclass setting on different folds in the case with 15 fully-annotated images and 213 weakly-annotated images.

		fold1	fold2	fold3	fold4	fold5	mean
	Standard supervision, whole tumor	74.31	78.91	67.57	75.55	75.96	74.46
	Mixed supervision, whole tumor	77.53 82.20	73.72 80.96 82.55 79.39
	Standard supervision, tumor core	61.17	63.89	55.72	55.36	63.18	59.87
	Mixed supervision, tumor core	62.83 65.26	62.23 61.99 67.23 63.91
	Standard supervision, enhancing core	66.15	64.83	53.83	61.68	62.77	61.85
	Mixed supervision, enhancing core	68.33 68.39	59.51 68.63 63.70 65.71

Preprint submitted to SPIE Journal of Medical ImagingDecember 10, 2018

Acknowledgements

Pawel Mlynarski is funded by the Microsoft Research-INRIA Joint Center, France. This work was supported by the Inria Sophia Antipolis -Méditerranée, "NEF" computation cluster.


The representation of Digital Factory Knowledge Entities by the Concept Map

Keywords: Digital Factory, Knowledge, knowledge management system "KSim", knowledge map

In a world characterized by increasing competitiveness, delivering a good product as soon as possible and increasing complexity of the product, the old method did not satisfy the needs of companies. In this context the Digital Factory (DF) is a solution that has emerged to design and simulate production systems throughout the product design process. In this case, different knowledge from different actors in different departments, are integrating in the Knowledge management system (KMS) named "KSim" and presented in the knowledge pages. With this variety of knowledge entities the actors spent a lot of time in the research of the relevant knowledge and the relationships between all knowledge entities to modify the files related to an error. We propose, in this paper, a module related to the KMS "KSim" able to represent all the knowledge entities used and obtained during the development of digital factory project in the conceptual map named "knowledge map".

Introduction

In the recent decades, a large number of companies sought to value the intangible investment (research and development, training, advertising, organizational methods, etc.) and in particular, their knowledge asset. This asset should be then re-used in different situations in order to reduce the costs and the time of development. Among such knowledge-based situations, Digital Factory (DF) has emerged to design and simulate production systems throughout the product design process. It can be defined as a set of software tools and methodologies allowing the design, simulation, initiation and optimization of production systems [1,[START_REF] Chryssolouris | Digital manufacturing: history,perspectives, and outlook[END_REF]. With this definition it is clear that knowledge management is very important in DF. In this context we developed a Knowledge management system (KMS) named "KSim" dedicated to the management of all knowledge entities in the Digital Factory [START_REF] Bouzid | KSim: an Information System for knowledge management in Digital Factory[END_REF]. In KSim, the knowledge is presented in knowledge pages that gather knowledge objects associated with a concept. Then we find in "KSim" a knowledge page of product that gathers all the products of a project, another knowledge page gathers all the parts of this project etc. however the actors found a great difficulty to see all the knowledge entities used and obtained during the Digital Factory project. With the complexity of product and the process of product development the actors spend a lot of time to research the relevant knowledge entities and to research the links between knowledge to correct the errors and modify all the files related to the error. For the purpose of solving all these problems, we develop in this paper a module related to "KSim" able to represent all the knowledge entities of Digital Factory project in the knowledge map to facilitate the understanding and the evolution of the project. The paper is structured as follows. Section 2 presents a study of the literature: first, we define knowledge entities then we analyse existing methods of knowledge representation. Section 3 describes the proposed framework. Section 4 describes the procedure of knowledge map construction with an example. Finally, Section 5 presents the conclusion.

State of the art

Knowledge Entities

Knowledge can be defined as refined, synthesized, systematized information, or as information associated with a context of use [START_REF] Baizet | Knowledge Management in design: Application to the computational mechanics at Renault-BED[END_REF][START_REF] Ayadi | Info Sim: Prototyping an information system for Digital Factory management[END_REF]. Merrill proposes a knowledge representation scheme consisting of knowledge components arranged into knowledge objects [START_REF] Merrill | Knowledge objects[END_REF]. The knowledge object or entity is a precise way to describe the subject matter content or knowledge to be taught. It is a framework for identifying necessary knowledge components. It is a way to organize a data base (knowledge base) of content resources (text, audio, video, and graphics) so that a given instructional algorithm (predesigned instructional strategy) can be used to teach a variety of different contents. Knowledge objects should consist of components that are not specific to a particular subject matter domain.

Knowledge representation method

In the literature, several techniques of knowledge representation have been devised for knowledge based applications in different domain insurance, banking, economic management, Process development Product…Among these methods we can find "Frame" [START_REF] Gordon | Creating knowledge maps by exploiting dependent relationships[END_REF], "Knowledge Page" that gathers knowledge entities associated with a concept, "Knowledge map" that a graphical representation of a set of knowledge entities… In this study, we are interested in the knowledge map which it is playing an important role in the KMS [START_REF] Yun | Knowledge-mapping model for construction project organizations[END_REF] to express the knowledge hierarchy and knowledge connection. Many research related to knowledge map have been carried in the recent decades. As an example [START_REF] Gómez | Knowledge maps: An essential technique for conceptualisation[END_REF] used a knowledge map to conceptualize the problem effectively and efficiently to gain a full understanding of a problem. [START_REF] Lv | Construction of intelligence knowledge map for complex product development[END_REF] proposed a knowledge map for the complex product design and development to improve the utilization efficiency of the knowledge using and shorten the time and effort spent on the Knowledge screening. [START_REF] Huang | Building a knowledge map model situated in product design[END_REF] proposed a knowledge map based on knowledge reuse process and its support to product design. But until now there is no common knowledge map able to represent all the knowledge used and developed during the development of Digital factory project and the link between these knowledge entities. So, we propose, in the next paragraph, to use a KMS named "KSim" [START_REF] Bouzid | KSim: an Information System for knowledge management in Digital Factory[END_REF], which able to manage all the corporate knowledge for the development of digital factory, to present the knowledge map of all knowledge entities used and developed in the digital factory project to facilitate the understanding of the project development process and to reduce the time of research of adequate knowledge entities.

Proposition

According to current conceptual model of KSim [START_REF] Bouzid | KSim: an Information System for knowledge management in Digital Factory[END_REF], the knowledge entities are linked. Indeed, the actors need to know and to see the general knowledge and need all the results obtained during this project to accomplish their works. Figure 1 shows the links between knowledge entities. For example the product is composed by several parts that have a materials and files of calculation, manufacturing and simulation by F.E. For the calculation file, the actor used a technical standard. Furthermore the experts of simulation need the results of product, process and resource to do the simulation…

Figure 1. The links between knowledge entities in KSim

To facilitate the understanding and evolution of the digital factory project and to know the relationships between the knowledge entities, we proposed to develop the knowledge map. To build it, we need three layers: graphic layer, description layer, and database layer (figure 2). The graphic layer refers the knowledge nodes and their relationships. The description layer is a link between the graphic layer and database layer. The database layer refers to the storage of all information and knowledge entities inserted in Info Sim [START_REF] Ayadi | Info Sim: Prototyping an information system for Digital Factory management[END_REF] and KSim. Knowledge map is made up knowledge nodes and knowledge links: • Knowledge node represents the name or title of knowledge entity and the name of creator or the person who has inserted this knowledge entity, • Knowledge link is the line among nodes represents the relationship between the nodes.

The procedure of knowledge map construction

To construct the knowledge map, we must first integrate the general knowledge in KSim and the project information in Info Sim [START_REF] Ayadi | Info Sim: Prototyping an information system for Digital Factory management[END_REF] and the relation between the knowledge entities. For example the figure 3 shows the technical standard that is a general knowledge entity used in the project knowledge of calculation. These integrated knowledge entities are stored in the database then presented in XML file that refers the description layer. The XML file is a source for the design of knowledge map. Finally, according to the data presented in figure 4 the knowledge map (figure 5) is made. Example of knowledge map Figure5 presents an example of knowledge map developed in our KMS "KSim" [START_REF] Bouzid | KSim: an Information System for knowledge management in Digital Factory[END_REF]. These figures show all the knowledge entities used and obtained during a project of Digital Factory and the relationships between all knowledge entities. For example in figure 5, we can see the links between simulation and resource, resource and process, process and product, product and part, part and files of calculation, calculation file and standard technical…

Conclusion

In this paper, starting from the definition of knowledge entity and the KMS "KSim", we developed a module for the design of knowledge map in "KSim". This module allows the presentation of all knowledge entities used and

Figure 2 .

 2 Figure 2. Conceptual layered model of knowledge map

Figure 3 .

 3 Figure 3. The knowledge description in KSim

Figure 4 .Figure 5 .

 45 Figure 4. Example of XML file

obtained during the development of Digital Factory and the relationships between them. The knowledge map helps the actors to understand and see the evolution of project, to know all the actors involved in this project, reduce the research time of relevant knowledge entity, see all the relationships between knowledge entities thus know the errors as soon as possible.

In the future work we will develop more this module to be able to present the characteristics of the links between knowledge (for example is composed, has…). Our KMS "KSim" able to support the integration of all actor viewpoints so we must think to develop a knowledge map able to represent all these viewpoints with their actors.

[1] W. Kuehn, "Digital Factory -Integration of simulation enhancing the product and production process towards operative control and optimisation," International Journal of Simulation, vol. 7, no. 7, p. 27-39, 2006.


Alexandre Gondran
email: alexandre.gondran@enac.fr

Laurent Moalic
email: laurent.moalic@uha.fr

Optimality Clue for Graph Coloring Problem

In this paper, we present a new approach which qualifies or not a solution found by a heuristic as a potential optimal solution. Our approach is based on the following observation: for a minimization problem, the number of admissible solutions decreases with the value of the objective function. For the Graph Coloring Problem (GCP), we confirm this observation and present a new way to prove optimality. This proof is based on the counting of the number of different k-colorings and the number of independent sets of a given graph G.

Exact solutions counting problems are difficult problems (#P-complete). However, we show that, using only randomized heuristics, it is possible to define an estimation of the upper bound of the number of kcolorings. This estimate has been calibrated on a large benchmark of graph instances for which the exact number of optimal k-colorings is known.

Our approach, called optimality clue, build a sample of k-colorings of a given graph by running many times one randomized heuristic on the same graph instance. We use the evolutionary algorithm HEAD [26], which is one of the most efficient heuristic for GCP.

Optimality clue matches with the standard definition of optimality on a wide number of instances of DIMACS and RBCII benchmarks where the optimality is known. Then, we show the clue of optimality for another set of graph instances.

Optimality Metaheuristics Near-optimal.

Introduction

For a given integer k ≥ 1, a k-coloring of a given graph G = (V, E) is an assignment of one of k distinct colors to each vertex v ∈ V in the graph, so that no two adjacent vertices (linked by an edge e ∈ E) are given the same color. The Graph Coloring Problem (GCP) is to find, for a given graph G, its chromatic number χ(G) corresponding to the smallest k such that there exists a k-coloring of G. GCP is NP-hard [START_REF] Karp | Reducibility among combinatorial problems[END_REF] for k ≥ 3. The k-coloring problem (k-CP) is the associated decision problem. For an optimization problem which is NP-hard, there is no efficient polynomial-time exact algorithm to solve it, unless P=NP. Therefore for large size instances of a minimization NP-hard problem, the exact algorithms must be stopped before their end. In this case, exact algorithms such as branch and bound methods find a lower bound of the optimal value of the objective function. Heuristic approaches are then the only ways to find, in reasonably fast running-time, a "good" solution in terms of objective function value, i.e. an upper bound of the optimal value. However, even if an admissible solution is found, its distance to the optimal solution remains unknown, except for approximation algorithms 1 . The optimality gap is the different between the upper bound (found by a heuristic) and the lower bound (found by a partial exact method). Optimality is proven only when this gap is equal to zero. Unfortunately for large size instances of an NP-hard problem, this gap is often important. It is particularly true for challenging instances [START_REF] Held | Maximum-weight stable sets and safe lower bounds for graph coloring[END_REF][START_REF] Moalic | Variations on memetic algorithms for graph coloring problems[END_REF] of the GCP of the DIMACS benchmark [START_REF]Cliques, Coloring, and Satisfiability: Second DIMACS Implementation Challenge[END_REF]. This paper addresses the following question: What to do in this situation? Is it possible to prove optimality of a graph coloring problem instance using only heuristic algorithms?

The response is Yes, for specific class of graphs: for example, it exists efficient polynomial-time exact algorithms to find χ(G) for interval graphs, chordal graphs, cographs [START_REF] Orlin | An O(n2) Algorithm for Coloring Proper Circular Arc Graphs[END_REF][START_REF] Shih | An O(n1.5) algorithm to color proper circular arcs[END_REF]. For some graphs like 1-perfect graphs 2 , for which the chromatic number χ(G) is equal to the size of the maximum cliques γ(G), it is possible to solve the dual problem, the Maximum Clique Problem (MCP), with another heuristic and conclude to optimality if the size of the maximum clique found is equal to the smallest number of colors used for coloring G found also by a heuristic. In this specific case, the optimality gap (or duality gap between GCP and MCP) is zero. However, the response to the question is No, in general case; a heuristic finds approximate solutions (upper bound); although the coloring found may be optimal, it is not possible to prove this possible optimality. Therefore, the question become: what can be done better using only a heuristic than finding an approximate solution? Is it possible to define a kind of optimality index for a graph coloring problem instance?

One shows in this article that a heuristic does not only find an upper bound of χ(G) but that it is also able to count the number of different k-colorings (i.e. the number of admissible solutions having the same objective function value). Our approach is based on the fact that the number of different k-colorings decreases dramatically when the number of colors, k, decreases too. Indeed figure 1 gives a typical example of a random graph with 30 vertices, a density of 0.9 and χ(G) = 16. The number of colorings with exactly k colors (blue bars) and the total colorings with k colors or less (red bars), noted N (G, k), are exactly computed for all values from k = 16 to k = 30. N (G, k) decreases exponentially when k decreases to χ(G). One proves a theorem showing that when the number of k-colorings is lower than a given value (the number of independent sets of G3 ), then we achieve the optimum: χ(G) = k.

In this article, we try to apply the proposed theorem in order to prove optimality.

Brief solutions counting review

Our work tackles the problem of counting solutions of NP-complete problems which has been widely studied for boolean SATisfiability problem, called #SAT, or Constraint Satisfaction Problem (CSP), called #CSP; kcoloring problem is a special case of CSP. These problems are known as #P-complete [START_REF] Valiant | The Complexity of Enumeration and Reliability Problems[END_REF]. A recent survey on #CSPs is done in [START_REF] Jerrum | Counting Constraint Satisfaction Problems[END_REF]. Even if a problem is not NP-hard, the problem of solutions counting is often hard. Specific studies on counting solutions of k-CP are done in [START_REF] Jerrum | A very simple algorithm for estimating the number of k-colorings of a low-degree graph[END_REF][START_REF] Favier | Solution counting for CSP and SAT with large tree-width[END_REF][START_REF] Miracle | Algorithms to approximately count and sample conforming colorings of graphs[END_REF]. Because the exact counting is in many cases a complex problem, statistical or approximate counting are often considered. Then, uniform sampling of the set of solutions problem is related to the problem of counting solutions. Many works are done on uniform or near uniform sampling like [START_REF] Gomes | From sampling to model counting[END_REF][START_REF] Gomes | Model Counting[END_REF][START_REF] Wei | A New Approach to Model Counting[END_REF]. The objective is to count by sampling. Frieze and Vigoda [START_REF] Frieze | A Survey on the use of Markov Chains to Randomly Sample Colourings[END_REF] give a survey on the use of Markov Chain Monte Carlo algorithms for approximately counting the number of k-colorings.

The features of ergodicity or quasi-ergodicity of the heuristics that guarantee an uniform sampling are deeply discussed in [START_REF] Ermon | Uniform Solution Sampling Using a Constraint Solver As an Oracle[END_REF]. However, theoretical results are obtained with a high value of k ≥ ∆ where ∆ is the maximum degree of the graph G which is very far from χ(G) for challenging graphs. On the other hand, when tests are performed with k = χ(G) like in [START_REF] Favier | Solution counting for CSP and SAT with large tree-width[END_REF], the considered graph instances are often with more than 10 20 k-colorings.

If the number of k-colorings is too high (higher than the number of independent sets), then it is not possible to apply our theorem. Therefore, in practice, our approach can be applied to graphs that do not have too many optimal colorings; we considered graphs with at most 1 million different optimal colorings. To our knowledge, it is the first time that solutions counting are used to prove optimality. We define a procedure, called optimality clue, in order to apply the proposed theorem. First, we build a sample of k-colorings of a given graph G by running many times (about 1,000 times) the same randomized heuristic algorithm. In this study, we use HEAD 4 , our open-source memetic algorithm (i.e. hybridization of tabu search and evolutionary algorithm), which is very efficient heuristic solving GCP [START_REF] Moalic | Variations on memetic algorithms for graph coloring problems[END_REF].

In this sample some colorings may appear several times and others only ones. The number of different kcolorings inside the sample is used to build an estimation of the total number of colorings with k colors. This estimator has been calibrated on a large benchmark of graph instances for which the number of optimal kcolorings is exactly known. Because we have no guarantee that the sampling is uniform, in the general case, therefore we have no guarantee that our estimator is always exact.

Moreover, building a sample of k-colorings is time-consuming, then the size of the sample should be "reasonable". Therefore, graphs for which our optimality clue can be calculated are graphs having not too many optimal k-colorings (i.e. about less than one million). Of course it is not possible to known a priori if a given graph has more or less than 1 million optimal colorings. Then, our approach provides a clue that a coloring found by the heuristic is perhaps optimal (positive conclusion) but never denies it (no negative conclusion): in many cases we can not have any conclusion.

This article is organized as follows. In Section 2, we present the new optimality proof for GCP based on solutions counting. Our general approach, called optimality clue, is define in Section 3. In Section 4, we detail how we calculate the estimate of the number of k-colorings using benchmark graph instances. Numerical tests and experiments are presented in Section 5. Finally, we conclude in Section 6.

Proof of Optimality by solutions counting

Notice that there are different ways to count the k-colorings of a given graph G. When counting the number of different k-colorings, we have to take into account the permutations of the color classes. We consider one k-coloring not as an assignment of one color among k to each vertex but as a partition of the vertices of the graph into k independent sets. An Independent Set (IS) or stable set is a set of vertices of G, no two of which are adjacent. Two k-colorings c1 and c2 are considered identical if they correspond to the same partition of G. The distance between two k-colorings that is taken into account is the set-theoretic partition distance used in [START_REF] Galinier | Hybrid evolutionary algorithms for graph coloring[END_REF][START_REF] Gusfield | Partition-distance: A problem and class of perfect graphs arising in clustering[END_REF][START_REF] Moalic | Variations on memetic algorithms for graph coloring problems[END_REF], which is independent of the permutation of the color classes. In previous works about solutions counting of k-CP [START_REF] Favier | Solution counting for CSP and SAT with large tree-width[END_REF], authors counted the total number of k-colorings including all the permutations like in the example of Figure 2; such a calculation of the number of different k-colorings is k! times higher than the way we count. This makes their methods inapplicable to our study. We write Ω(G, k) the set of all k-colorings of the graph G. A k-coloring can use exactly k colors or less, then

Ω(G, k -1) ⊂ Ω(G, k). The cardinal of Ω(G, k) is noted N (G, k) = |Ω(G, k)|.
Our approach is based on the following fact :

Lemma 1 Let a graph G and an integer k ≥ 1. If there exists at least one k-coloring of G, then there exists at least i(G) -k + 1 different (k + 1)-colorings of G:

N (G, k + 1) ≥ i(G) -k + 1,
where i(G) is the number of independent sets of G.

Proof Notice that a k-coloring of a graph G = (V, E) is a partition of |V | vertices into k IS. Indeed vertices colored with the same color inside a k-coloring are necessarily an IS. In other words, it is always possible to color all vertices of any IS with the same color. We note

IS(G) = {U ⊂ V | ∀x, y ∈ U 2 , {x, y} / ∈ E} the set of all the IS of G, then i(G) = |IS(G)|.
Starting with one coloring of G with exactly k colors, for each independent set of G except for the k IS of the k-coloring, it is possible to recolor all vertices of this independent set with a new color (the (k + 1)th color). We obtain by this way one different (k + 1)-coloring for each different independent set, then we count at least a total of i(G) -k different colorings with exactly (k + 1) colors. Then, N (G, k + 1) ≥ i(G) -k + 1 because we have to count also the starting k-coloring.

Then, we obtain the following theorem:

Theorem 1 Let a graph G and an integer k ≥ 1. Let N (G, k) the number of k-colorings of G and i(G) the number of independent sets of G. If i(G) -k > N (G, k) > 0, then χ(G) = k. Proof χ(G) ≤ k because N (G, k) > 0. If χ(G) < k
, it means that there exists at least one (k -1)-coloring (i.e. N (G, k -1) > 0). If we add a new color, it is possible to consider this (k -1)-coloring and to recolor any independent set of G with the new color, we obtain by this way i(G) -k different k-colorings (by Lemma 1). Therefore i(G) -k ≤ N (G, k) which refute initial assumption.

For example, the studied graph in Figure 1 (30 vertices and density 0.9) has 38 different colorings with 16 colors: N (G, k = 16) = 38; moreover this graph has 78 IS: i(G) = 78, then the theorem is applicable with k = 16 because: i(G) -k = 78 -16 = 62 > 38 = N (G, k) > 0. Then, thanks to the theorem we can conclude that χ(G) = 16. Moreover, for k = 17, N (G, k = 17) = 3121 > i(G) -k = 61, so the theorem is not applicable.

Corollary 1 Let a graph G and an integer

k ≥ 1. Let N (G, k) an upper bound of the number of N (G, k) and i(G) a lower bound of i(G). If N (G, k) > 0 and i(G) -k > N (G, k), then χ(G) = k.

Optimality Clue

We propose in this paper to apply the corollary 1, so to find an appropriate upper bound of the number of k-colorings of G, N (G, k), and a lower bound of the number of independent sets of G, i(G).

IS counting

There exists many algorithms [START_REF] Bron | Algorithm 457: Finding All Cliques of an Undirected Graph[END_REF][START_REF] Carraghan | An exact algorithm for the maximum clique problem[END_REF][START_REF] Östergård | A fast algorithm for the maximum clique problem[END_REF][START_REF] Samotij | Counting independent sets in graphs[END_REF] for counting all the maximal independent sets of a graph G (or similarly counting all the maximal cliques5 in G, the complementary graph of G). By definition, the number of maximal IS, noted imax(G), is a lower bound of i(G). Those algorithms are based on enumeration. Because we focus this study on graphs having less than 1 million optimal solutions, we can stop the enumerating after finding 1 million IS. Generally, i(G) is very high except for graphs with very high density. Real-life graphs have often a low density, then i(G) is very high. Moreover, a simple lower bound is given by [START_REF] Pedersen | Bounds on the Number of Vertex Independent Sets in a Graph[END_REF] :

i(G) ≥ 2 α(G) + n -α(G),
where α(G) is the size of the largest independent set of G and n the number of vertices. Bollobás' book [START_REF] Bollobás | Random Graphs[END_REF] (p.283) gives also a statistical number of maximal cliques of size p for a random graph. Then, we conclude that:

imax(G) ≈ iB(G) = n p=1 n p (1 -d) ( p 2 )
with n the number of vertices and d the density of a random graph G.

In this study, we use Cliquer6 , an exact branch-and-bound algorithm developed by Patric Östergård [START_REF] Östergård | A fast algorithm for the maximum clique problem[END_REF] that enumerates all cliques (an IS is a clique in the complementary graph).

It is more complex to evaluate N (G, k) and section 4 presents a way to build an experimental upper bound of N (G, k). We characterize this upper bound as experimental because it is based on experimental tests on benchmark graph instances, then there is no total guaranty that it is an upper bound.

Procedure

We define here the procedure of what we call Optimality Clue for graph coloring: let G a graph and k > 0 a positive integer, that we suspect to be the chromatic number of G. The proposed approach is based on the five following steps:

1. Build a sample of t = 1, 000 k-colorings of G: we run the memetic algorithm HEAD on G as many times as needed to obtain t legal k-colorings. Those solutions are the solutions sample. The size of the sample is equal to t. We take in general case t = 1, 000 when it is possible.

(G, k) = 100000 (G, k) = 500000 (G, k) = 1000000
Figure 4: Collision probability q, given the sample size t, and the total number of k-colorings N (G, k).

Uniform sample

If the sample is uniform7 , then there exists statistical methods to count solutions and to build an upper bound with statistical guarantee, for example the capture-recapture methods: Peterson method [START_REF] Krebs | Ecology[END_REF], Jolly-Seber method [START_REF] Baillargeon | Rcapture: Loglinear Models for Capture-Recapture in R[END_REF] wich is commonly used in ecology to estimate an animal population's size. However, it is not our case: we have no guarantee that our solutions sample is uniform or near uniform. HEAD is a memetic algorithm that explores the space of non-legal k-colorings: a non-legal k-coloring is a coloring with at most k colors and where two adjacent vertices (linked by an edge) may have the same color (called conflicting edge).

The objective of HEAD is to minimize the number of conflicting edges to zero, that is to get a legal k-coloring.

HEAD is an evolutionary algorithm with a population size equals to two. The two non-legal k-colorings perform at each generation a tabu search and after a crossover. The sample distribution depends on the fitness landscape properties [START_REF] Merz | Memetic algorithms for combinatorial optimization problems: fitness landscapes and effective search strategies[END_REF][START_REF] Marmion | Fitness Landscape Analysis and Metaheuristics Efficiency[END_REF] 8 and there is no reason for this distribution to be uniform. A smooth landscape (respectively a rugged landscape) around a legal k-coloring will increase (resp. decrease) the probability of finding this k-coloring. Figure 3 represents the frequency of the 319 optimal 46-colorings of <r140_90.4> graph of RCBII benchmark (140 vertices and density 0.9) in a sample of size 100,000 found by HEAD heuristic. In this typical graph instance, the ratio between the least frequent and the most frequently found coloring is around a factor of 10 3 which corresponds to the same scale as similar studies [START_REF] Wei | A New Approach to Model Counting[END_REF].

Another approach is to take into account the ergodicity of an algorithm, which is its capability to explore all the search space. More precisely, an algorithm is ergodic if it is possible (probability not null) to reach any k-coloring from any other k-coloring in a finite number of iterations. Random walks or Metropolis algorithms (with a positive temperature sufficiently high) are ergodic algorithms since there is always a finite probability of escaping from local minimum. However, those algorithms are very inefficient in practice to find an optimal k-coloring in the general case.

Sample size

The choice of t, the size of the sample, is very important for two reasons. First, in practice, to build a sample of k-colorings can be very time-consuming, then the size of the sample should have a reasonable size. We take t = 1, 000 for most of the graph instances. However, the more challenging the graph instance, the longer HEAD takes to find one k-coloring. Therefore, it is not possible to build a sample of size 1,000 for all graphs, such as for the <DSJC500.5> graph of DIMACS (cf. Table 2).

The second reason is more theoretical. We have limited the maximum number of different optimal solutions to 1 million, for a graph to be considered by our approach. In fact, we choose 1 million because it equals to t 2 with t = 1, 000. Indeed, if the sample is uniformly drawn at random in Ω(G, k), the probability q that at least two colorings of the sample are identical is equal to 9

: q = 1 - N ! N t (N -t)! 1 -e -t(t-1)

2N

then N -t(t-1) 2ln(1-q) . We call also q the collision probability. So, if q = 0.5 then N 720626, if q = 0.393 then N ∼ t 2 = 10 6 . Figure 4 represents the collision frequency, q, in function of the sample size, t, for different values of the Ω(G, k) size.

When N (G, k) = 10 5 and t = 1, 000, it is almost impossible to miss a collision in the sample, but for N = 10 6 , there is around 60% to miss a collision. However, it is not tragic to miss a collision for our approach. Indeed, the consequence is that the clue of optimality may be not applicable but the risk of false positive is avoided. A false positive occurs if our procedure 3.2 improperly indicates the optimality clue, when in reality the k-colorings are not optimal. Moreover, the collision frequency is higher for a non-uniform sample than for a uniform one.

4 Estimate of the number of k-colorings: U B(G, k, p, t)

Data sets

In order to define an estimator or at least an upper bound of the number of k-colorings, we need to have a large number of graph instances for which we know the exact number of k-colorings. Fabio Furini et al. [START_REF] Furini | An Improved DSATUR-Based Branch-and-Bound Algorithm for the Vertex Coloring Problem[END_REF] have published an open-source and very efficient version of the backtracking DSATUR algorithm [START_REF] Brélaz | New Methods to Color the Vertices of a Graph[END_REF] which returns the chromatic number of a given graph 10 . DSATUR is one of the best exact algorithms for GCP, particularly for graphs with high density. We suggest readers interrested in an overview of exact methods for GCP to read [START_REF] Malaguti | A survey on vertex coloring problems[END_REF][START_REF] Held | Maximum-weight stable sets and safe lower bounds for graph coloring[END_REF].

We modified their DSATUR algorithm in order to count the total number of k-colorings. The pseudo code of the algorithm, called CDSATUR, is presented in algorithm 1. CDSATUR returns, for all values k, the exact value of N (G, k) taking into account the permutation of colors and especially N (G, k = χ(G)).

Algorithm 1 CDSATUR which returns the number of all k-colorings of G: N (G, k). Data: G = (V, E) a graph and k a positive integer.

N ← 0 C[v] ← N one, ∀v ∈ V : C is the empty coloring. l ← 0: number of colors used by C. CDSAT U R(C, l) return N Procedure CDSAT U R(C , k ): if all the vertices of C are colored then if k ≤ k then N ← N + 1 else Select an uncolored vertex v of C for every feasible color i ∈ [1 ; k + 1] do C ← C , C [v] ← i k ← max(k , i): number of colors used in C . if k ≤ k then CDSAT U R(C , k )
Fabio Furini et al. published also 2031 random GCP instances called RCBII 11 with vertices from 60 to 140 and density between 0.1 and 0.9. This wide variety of graphs is our reference dataset. We complete this dataset with easy DIMACS graphs [START_REF]Cliques, Coloring, and Satisfiability: Second DIMACS Implementation Challenge[END_REF] for which χ(G) and N (G, χ) is computable with CDSATUR.

The 2031 graphs of RCBII benchmark have characteristics described in Table 1. We can notice that χ(G) is known for all these graphs [START_REF] Furini | An Improved DSATUR-Based Branch-and-Bound Algorithm for the Vertex Coloring Problem[END_REF]. First we calculated N (G, χ) with CDSATUR, with a time limit equals to 2400s. This time is enough for most of the graphs. There are only 210 graph instances of RBCII (on the 2031) for which CDSATUR does not have enough time to find N (G, χ). These 210 graphs are used to test our approach (test dataset).

Among the graphs for which N (G, χ) can be determined, we consider only those with less than 1 million optimal solutions: they form the reference dataset (959 graph instances). Finally, we can distinct inside the reference dataset, graph instances verifying i(G) > N (G, χ) (566/959) or not (393/959).

It remains 862 graphs on the 2031 of RBCII benchmark with more than 1 million of optimal solutions. We decided to test our approach on those graphs (called control dataset) to check if the proposed algorithm can produce false positives or not.

Analysis of graph instances

Before determining an upper bound of N (G, χ), we investigate the possible links between standard features of a graph as its size (number of vertices), its density, or its chromatic number and the number of optimal colorings: N (G, χ)

Links between N (G, χ), graph size and density

Graphs with same size (number of vertices) and same density can have a number of optimal colorings very different from one another. A typical example is given in Figure 5 where is represented the distribution of 49 graph instances with 80 vertices and density 0.3 (<r80_30.*> of RBCII benchmark) in function of the number of solutions N (G, χ). Half of the graphs (25/49) have less than 100 000 optimal solutions while a third (18/49) have more than 1 million optimal solutions. There are no simple law that characterize this distribution. However, we can notice that the lower the density, the higher the optimal solution number. Indeed, Figure 6 presents the proportion of graphs with 70 vertices of RBCII benchmark having more than 1 million colorings depending on graph density. For a low density such as 0.1, nearly all graphs have more than 1 million optimal solutions, while no graph with high density (equals to 0.9).

In order to have a more fine view of the link between the number of optimal colorings and the graph density, we generated 1,000 random graphs with 50 vertices and density d (d = 0.1, 0.2, ..., or 0.9). Each line in Figure 7 represents (for each density) the proportion of graphs having less than n optimal colorings with n between 10 2 and 10 6 . Pink line of Figure 7 shows for example that 50% of graphs (with 50 vertices and density = 0.3) have less than 10 5 optimal colorings. The plots are quite similar for graphs with 60 or 70 vertices. The graph size seams to have a slight influence on the number of optimal solutions.

Links between N (G, χ) and χ(G)

As shown in Figure 8, there is no obvious link between the chromatic number of a graph, χ(G) (y-axis) and the number of optimal colorings N (G, χ) (x-axis). Each dot of Figure 8 corresponds to one graph of RCBII for which it is possible to calculate exactly N (G, χ) with CDSATUR. Figure 9: Each blue dot corresponds to one of the 959 graph instances that have less than 1 million optimal colorings (reference dataset). The number of optimal colorings (calculated with CDSATUR algorithm) is in ordinate. The number p of different optimal solutions found by our HEAD algorithm after 1,000 success runs is in abscissa. The red line is an upper bound U B(G, k, p, t) of the number of optimal colorings. The right figure is a zoom of left figure for p ≤ 500.

Upper bound function

of different colorings p found by HEAD on the total of t = 1, 000 success runs (in abscissa) and the exact number of colorings, N (G, k), calculated with CDSATUR (in ordinate). Each dot corresponds to one graph of the reference dataset. The objective now is to determine an upper bound of N (G, k), U B, as small as possible. Indeed, in order to apply the Theorem 1, we must have i(G) -k > U B. Figure 9-right which is a zoom of the left figure for p ≤ 500 shows that for p t, p is near linear to N (G, k): p ∼ N (G, k). Then, p is a good candidate to be an estimator of N (G, k). When p is near to t, the range of N (G, k) values is very large, near to p 2 = 10 6 , and p is a very bad estimation of N (G, k), but notice that N (G, k) < p 2 . We add on those figures a red line that represents a possible upper bound of N (G, k) which is equal to:

U B(G, k, p, t) = p + p α t+p t if p < t × 0.99 +∞ otherwise (1)
with α = 1.01. Indeed, when p t, U B(G, k, p, t) ∼ 2p and when p is near to t, U B(G, k, p, t) ∼ p 2 . Between these extreme values, the cloud of blue dots follows approximately an exponential curve. U B(G, k, p, t) was also built to be above all blue dots; i.e. it is a valid upper bound for all graphs of the reference dataset. Of course, there is no guarantee that this upper bound is still valid for all other graphs. So, our approach is never able to prove optimality in a strict sense. It gives only a clue.

Experiments and analysis

Tests

The upper bound U B was built based on the graphs of the reference dataset. Now, in order to test the optimality clue (procedure Section 3.2), we use this upper bound on graphs of the test dataset and the control dataset and for some graphs coming from the DIMACS benchmark.

Results on RCBII benchmark are presented in Table 1 in the two last lines. The first column concerns the 862 graphs with more than 1 million optimal solutions, corresponding to the control dataset. There is no false positive: the procedure 3.2 concludes for all the graphs that there is no optimality clue. The two following columns concern the reference dataset. More precisely, the second column concerns graphs having less than 1 million optimal solutions but that do not verify Theorem 1: the number of IS is lower than the number of optimal solutions. Of course, there are no false positives for this case, because U B was built to validate those graphs (reference dataset). The third column concerns the 566 graphs verifying the Theorem 1. The optimality clue is proven for 449 of them because i(G) > U B(G, k, p, t) > N (G, k). The optimality clue is not shown on the 117 (= 566 -449) other graphs because U B(G, k, p, t) ≥ i(G) > N (G, k). U B(G, k, p, t) is an upper bound too high in this case. To prove the optimality clue on those graphs, we would have to increase the size of the solutions sample, t. The fourth column concerns the test dataset i.e. graphs for which the number of optimal solutions is unknown. We prove the optimality clue for nearly 20% of these graphs (39/210). There are three reasons why we did not prove the optimality clue for the other 171 (=210-39) graphs: 1) graph instances have more than 1 million solutions; 2) graph instances do not verify the Theorem 1; Nothing can be done for these two first reasons. 3) p is too close from t, then the upper bound U B is too high. In order to have an upper bound more accurate, i.e. still valid but not too high, we have to increase the size of the sample or to choose another formula than equation [START_REF] Baillargeon | Rcapture: Loglinear Models for Capture-Recapture in R[END_REF]. Our approach therefore applies to about 20% of the random graphs in the RCBII benchmark. For control and reference datasets, we get more or less the same proportion: 25% (449/1821).

Table 2: Results of optimality clue tests for graphs of DIMACS benchmark with p < t. The results on selected DIMACS benchmark graphs are presented in Table 2. We only present graph instances for which the solutions sample generated by HEAD, are not all different (i.e. p < t = 1, 000) and are susceptible to be optimal.

Instances |V | d χ(G) k i(G) N (G, k) t p U B(G,
The first column of Table 2 indicates the name of the graph instance. Columns 2-7 indicate for each graph, its number of vertices |V |, its density d, its chromatic number χ(G), when it is known, the number of colors k, used for the test (k = χ(G) if χ(G) is known), its number of independent sets i(G), and the exact number of legal k-colorings N (G, k), when it is possible to calculate it with CDSATUR. Then, columns 8-10 indicate the size of the solutions sample t, the number of different solutions in the sample p, the experimental upper bound of the number of k-colorings U B(G, k, t, p). Columns 11 and 12 provide the result of the procedure of Section 3.2 of optimality clue and the total computation time in seconds to generate all the sample. The two last columns indicates the lower bound of χ(G) found by the best known exact method [START_REF] Held | Maximum-weight stable sets and safe lower bounds for graph coloring[END_REF] (or by IncMaxCLQ [START_REF] Li | Combining maxsat reasoning and incremental upper bound for the maximum clique problem[END_REF], that found the maximum clique) and the computation time of this method.

The first part of Table 2 corresponds to 17 graphs for which χ(G) is already known by other methods. We prove the optimality clue for 12 of them. Computation time of optimality clue is higher than those for finding the upper bound with the exact methods except for two graphs. However the computation time of optimality clue can be considerably reduced because the 1,000 runs of HEAD can be switch on 2, 3... or 1,000 different processors or on a cluster of computers.

The second part of Table 2 (below the horizontal line) corresponds to 6 graphs for which χ(G) is unknown. We have generated 3 new graph instances called <DSJC+*.*_k> with rules almost similar of those of <DSJC*.*> but for which a k-coloring is hidden and so that very few others k-colorings can exist. For <DSJC500.5> which is one of the most challenging graph of DIMACS, we have the optimality clue for k = 47 12 .

Notice that when the optimality clue is proven for a given k, we check that the optimality clue for k + 1 can not be proved as well. For example, for <DSJC500.5>, we check the optimality clue for k = 48 is not achieved: for t = 100, 000 48-colorings found by HEAD, all are different. If the randomized algorithm HEAD is biased, i.e. finds k-colorings always in the same subset of Ω(G, k) (and therefore undervalues the upper bound), this bias does not reveal for k + 1.

Notice that the impact of the size of the sample has a great importance on the test. For having a not too high upper bound of N (G, k), we should have p t. So, the size t of the solutions sample can be chosen in function of the results of each graph instance. For example, we can extend the sample until all colorings of the sample are found at least twice (cf. Good-Turing estimator).

Conclusions and perspectives

Based on Theorem 1, we propose a procedure, called optimality clue, for determining if the global optimum is reached or not, by a heuristic method.This approach estimates an upper bound of the number of legal k-colorings by running a randomized heuristic several times. This process is contextual to the instance to solve. No general conclusion can be drawn on the heuristic itself, which is used to build solutions. This definition can be seen as an experimental criterion that evaluates the convergence of a randomized algorithm to the chromatic number. However, since it is not possible to be sure that the upper bound is exact, it is not possible to prove optimality in the strict sense.

Our approach is nevertheless an alternative when the exact methods are not applicable (high optimality gap). It is a new way for providing a criterion on the proximity of the optimality. The general idea is that the number of solutions with a same objective function value decrease when the objective function is getting closer to the optimal value. Optimality clue matches with the standard definition of optimality on a wide number of instances of DIMACS and RBCII benchmarks where the optimality is known. Furthermore, we proved the optimality clue for <DSJC500.5> graph of DIMACS with k = 47 colors which is a very challenging instance (only two algorithms are able to find 47-colorings [START_REF] Titiloye | Parameter Tuning Patterns for Random Graph Coloring with Quantum Annealing[END_REF][START_REF] Moalic | Variations on memetic algorithms for graph coloring problems[END_REF]). Tests on small random graphs (under 140 vertices) show that optimality clue can be proved for 20% of them.

Finally, we defined an upper bound quite high to avoid false positives: graphs for which we prove the optimality clue for a given k, while χ(G) = k. 13 .

Representative sample

The proposed approach is based on a sampling of the legal k-colorings space, Ω(G, k). This sampling is built by running many times HEAD algorithm; each success run providing one element of the sample. Ideally, to obtain a representative sample, HEAD has to uniformly draw one k-coloring inside the legal k-colorings space. Of course, it is not possible to guarantee this feature in all cases, it is why we built an upper bound function (Equation 1) of N (G, k) = |Ω(G, k)|. In order to improve our approach and to get closer to the ergodic objective, we plan to use more powerful model counting such as presented in [START_REF] Gomes | Model Counting[END_REF][START_REF] Ermon | Uniform Solution Sampling Using a Constraint Solver As an Oracle[END_REF] and study the ergodic propriety of HEAD. Our work is only an initial contribution to the study of optimality by counting. Other methods of estimating the population should be tested, such as Good-Turing methods that estimate missing mass (i.e. missing k-colorings in the sample) or Peterson-type methods to obtain statistical guarantees.

Generalization to other optimization problems All tests presented in this paper are done on GCP, which has the special propriety of Lemma 1. Then, it is possible to generalize our approach to other problems, as soon as they have an analogue propriety of Lemma 1. This is the case for the maximum clique problem.

Appendix: Generalization to other optimization problems

Here is a presentation of the conditions to be fulfilled for a generalization of the proposed approach to other optimization problems. Then, we recall some basic notations. An optimization problem can always be modeled as follows:

f, X Minimize f (x) s.c. x ∈ X
where X is the set of admissible solutions, x the decision variable and f : X → R is the objective function. We note f, X this problem, S( f, X , k) = {x ∈ X|f (x) = k} the set of admissible solutions with an objective function value equal to k and N ( f, X , k) = |S( f, X , k)| the size of this set, i.e. the number of admissible solutions with objective function value equal to k. For the GCP, X corresponds to the set of legal colorings and f provides the number of colors used.

Definition 1 Knowing one instance of an optimization problem (minimization case) f, X and k a real number. If:

1. (x1, ...xt) are t admissible solutions drawn independently from S( f, X , k); it is a sample of S( f, X , k),

2.
N is an estimator of the upper bound of N ( f, X , k), the number of different solutions with value equal to k;this estimator is based on the sample (x1, ...xt).

3. lb > 0 is a positive integer such as if ∃ε > 0, S( f, X , k -ε) = ∅ then N ( f, X , k) ≥ lb; that is to say, if there exists at least one solution of the optimization problem f, X with better objective function than k, then the number of solutions of S( f, X , k) is higher than lb. This is analogous to Lemma 1 for GCP.

lb > N

then we say that solutions of S( f, X , k) have a optimality clue of the optimization problem f, X relative to the upper bound N .

To estimate N , it is possible to adopt the same method as in Section 5. However, defining lb can be quite complex or impossible depending of the problem. For example, for the Traveling Salesman Problem (TSP), if it exits a solution with a length equal to n -then, how is it possible to estimate the number of solutions with length equals to n?

Optimaly clue can be easyly tested for example on Maximum Clique Problem (MCP). If we suspects that the maximum clique of a given graph is γ(G) = 29. A simple lower bound can be determinated. Indeed, if a graph has an unique maximum clique of size 30, for example, then there exists at least 29 cliques with size 29, i.e. : lb = γ(G) = 29. This value may be small but the optimality is proved if lb > N with N the total number of clique of size 29.

Figure 1 :

 1 Figure 1: Number of colorings with exactly k colors (blue bars) and number of total colorings with k colors or less (red bars), noted N (G, k) in function of k = 16...30, for a random graph with 30 vertices, density 0.9 and χ(G) = 16.

Figure 2 :

 2 Figure 2: Two 3-colorings c1 and c2 of the same graph with four vertices. These two colorings have to be considered identical because d(c1, c2) = 0 with d the set-theoretic partition distance; we pass from one to the other, just by a permutation of color classes.

Figure 5 :

 5 Figure 5: Histogram characterizing the distribution of the 49 random graphs with 80 vertices and a density of 0.3 (<r80_30.*> for the RBCII benchmark) given the number of optimal colorings, N (G, χ).

  of graph <r_70.*.*> to have more than 1 millions of optimal colorings

Figure 6 :

 6 Figure 6: Proportion of graphs with 70 vertices (<r70_*.*> of RBCII benchmark) having more than 1 million of optimal colorings given the density.

  We define in this Section an upper bound of N (G, k) based on the 953 graphs of the reference dataset. Suppose we have, for a given graph G, a set of n different k-colorings:Ω(G, k) = {x1, ..., xn}, i.e. n = |Ω(G, k)| = N (G, k)is unknown. We also have a sequence W of t independent samples: W = (w1, ..., wt), where w k ∈ Ω(G, k), ∀k = 1...t. This sample W is composed of t independent success runs of HEAD algorithm. We note ∀j = 1...n, #(xj) the count of xj in W . For these t colorings, we count p different colorings in W : p = |{xj ∈ W, #(xj) > 0}|. So then, N (G, k) ≥ p and t ≥ p ≥ 1. Figures 9 represent for each graph of the reference dataset, the number

Figure 7 :Figure 8 :

 78 Figure 7: Proportion of random graphs (with 50 vertices and a given density) having less than n optimal colorings with n range between 1 hundred to 1 million.

Table 1 :

 1 Distribution of 2031 RCBII graph instances

							χ known	Total
			N > 10 6			N ≤ 10 6	N ?	#instances
	862 (control dataset)		959 (reference dataset)	210 (test dataset)	2031
						i(G) ≤ N	i(G) > N
							393		566
	opt. clue	not opt. clue	opt. clue	not opt. clue	opt. clue	not opt. clue	opt. clue	not opt. clue
		0		862		0	393	449	117	39	171
		18
		16
		14
		12
	Nb of graphs	8 10
		6
		4
		2
		0	0	200 000	400 000 Nb of optimal solutions 600 000	800 000 < 1 000 000

Notice that it is still NP-hard to approximate χ(G) within n 1-for any > 0[START_REF] Zuckerman | Linear degree extractors and the inapproximability of max clique and chromatic number[END_REF].

A perfect graph is a graph in which the chromatic number of every induced subgraph equals the size of the largest clique of that subgraph. 1-perfect graphs are more general than perfect graphs. There exists polynomial-time exact algorithms to find χ(G) for perfect graphs[START_REF] Grötschel | Polynomial Algorithms for Perfect Graphs[END_REF], but slow in practice. Line graphs, chordal graphs, interval graphs or cographs are subclasses of perfect graphs.

An independent set is a subset of vertices of G, such that every two distinct vertices in the independent set are no adjacent.

Open-source code available at: github.com/graphcoloring/HEAD

A maximal clique is a clique that cannot be extended by including one more adjacent vertex. A maximum clique is a clique that has the largest size in a given graph; a maximum clique is therefore always maximal, but the converse does not hold. Analogue definition for IS.

Code available at: users.aalto.fi/ pat/cliquer.html. To count all IS of a graph, you just execute: ./cl <complement graph> -a -m 1 -M <k>

All k-colorings in the sample are uniformly drawn at random in Ω(G, k).

The fitness landscape itself depends on the neighborhood used for tabu search and the crossover used.

This problem is linked to the birthday problem that shows that in a room of just 23 people there's a 50-50 chance that two people have the same birthday. In our case, the number of days in a year is N and the number of people is the size t of the sample.

Code available at: lamsade.dauphine.fr/coloring/doku.php

Instances available in the same address

For <DSJC500.5> the computation time is not report because it takes several weeks and no accurate time has been recorded.

In this context, we propose on our website a challenge to find a counterexample (false positive graph)