Hao Hu

Marie-José Huguet
email: huguet@laas.fr

Mohamed Siala
email: siala@laas.fr

Optimizing Binary Decision Diagrams with MaxSAT for Classification

The growing interest in explainable artificial intelligence (XAI) for critical decision making motivates the need for interpretable machine learning (ML) models. In fact, due to their structure (especially with small sizes), these models are inherently understandable by humans. Recently, several exact methods for computing such models are proposed to overcome weaknesses of traditional heuristic methods by providing more compact models or better prediction quality. Despite their compressed representation of Boolean functions, Binary decision diagrams (BDDs) did not gain enough interest as other interpretable ML models. In this paper, we first propose SAT-based models for learning optimal BDDs (in terms of the number of features) that classify all input examples. Then, we lift the encoding to a MaxSAT model to learn optimal BDDs in limited depths, that maximize the number of examples correctly classified. Finally, we tackle the fragmentation problem by introducing a method to merge compatible subtrees for the BDDs found via the MaxSAT model. Our empirical study shows clear benefits of the proposed approach in terms of prediction quality and interpretability (i.e., lighter size) compared to the state-of-the-art approaches.

Introduction

Due to the increasing concerns in understanding the reasoning behind AI decisions for critical applications, interpretable Machine Learning (ML) models gained a lot of attention. Examples of such ML applications include job recruitment, bank credit applications, and justice [START_REF] Voigt | The EU General Data Protection Regulation (GDPR): A Practical Guide[END_REF]. Most of traditional approaches for building interpretable models are greedy, for example, decision trees [START_REF] Breiman | Classification and Regression Trees[END_REF][START_REF] Quinlan | Induction of Decision Trees[END_REF][START_REF] Quinlan | C4.5: programs for machine learning[END_REF], rule lists [START_REF] Cohen | Fast Effective Rule Induction[END_REF][START_REF] Clark | Rule induction with CN2: Some recent improvements[END_REF], and decision sets [START_REF] Kaufmann | Interpretable Decision Sets: A Joint Framework for Description and Prediction[END_REF]. Compared to traditional approaches, exact methods offer guarantee of optimality, such as model size and accuracy. In this context, combinatorial optimisation methods, such as Constraint Programming [START_REF] Bonfietti | Embedding Decision Trees and Random Forests in Constraint Programming[END_REF][START_REF] Verhaeghe | Learning optimal decision trees using constraint programming[END_REF], Mixed Integer Programming [START_REF] Angelino | Learning certifiably optimal rule lists for categorical data[END_REF][START_REF] Verwer | Learning optimal classification trees using a binary linear program formulation[END_REF][START_REF] Aglin | Learning Optimal Decision Trees Using Caching Branch-and-Bound Search[END_REF], or Boolean Satisfiablility (SAT) [START_REF] Bessiere | Minimising Decision Tree Size as Combinatorial Optimisation[END_REF][START_REF] Narodytska | Learning Optimal Decision Trees with SAT[END_REF][START_REF] Avellaneda | Efficient Inference of Optimal Decision Trees[END_REF][START_REF] Hu | Learning Optimal Decision Trees with MaxSAT and its Integration in AdaBoost[END_REF][START_REF] Janota | SAT-Based Encodings for Optimal Decision Trees with Explicit Paths[END_REF][START_REF] Yu | Computing Optimal Decision Sets with SAT[END_REF]) have been successfully used to learn interpretable models. These declarative approaches are particularly interesting since they offer certain flexibility to handle additional requirements when learning a model.

By providing compact representations for Boolean functions, Binary Decision Diagrams (BDDs) [START_REF] Akers | Binary Decision Diagrams[END_REF][START_REF] Moret | Decision Trees and Diagrams[END_REF][START_REF] Bryant | Graph-Based Algorithms for Boolean Function Manipulation[END_REF][START_REF] Knuth | The Art of Computer Programming, Volume 4, Fascicle 1: Bitwise Tricks & Techniques; Binary Decision Diagrams[END_REF]) are widely studied for hardware design, model checking, and knowledge representation. In the context of ML, BDD could be viewed as an intrepretable model for binary classification. In addition, they were extended for multi-classification, known as decision graphs and heuristic methods were proposed in [START_REF] Oliver | Decision graphs: an extension of decision trees[END_REF][START_REF] Kohavi | Bottom-Up Induction of Oblivious Read-Once Decision Graphs[END_REF][START_REF] Kohavi | Oblivious Decision Trees, Graphs, and Top-Down Pruning[END_REF][START_REF] Mues | Decision Diagrams in Machine Learning: An Empirical Study on Real-Life Credit-Risk Data[END_REF]). Moreover, [START_REF] Ignatov | Decision Stream: Cultivating Deep Decision Trees[END_REF] proposed decision stream, a similar topology to BDD based on merging similar subtrees in each split made in decision trees to improve the generalization. [START_REF] Oliver | Decision graphs: an extension of decision trees[END_REF][START_REF] Kohavi | Bottom-Up Induction of Oblivious Read-Once Decision Graphs[END_REF] showed that decision graphs could avoid the replication problem and fragmentation problem of decision trees effectively, which BDDs also could avoid in binary classification. This fact indicates that generally in the practice of ML, a BDD have a smaller size than the corresponding decision tree.

In this paper, we introduce a SAT-based model for learning optimal BDDs with the smallest number of features classifying all examples correctly, and a lifted MaxSAT-based model to learn optimal BDDs minimizing the classification error. We assume that all BDDs are ordered and reduced 1 , the limitation on the depth for a BDD, corresponds to the number of features to be selected by our model. To the best of our knowledge, [START_REF] Cabodi | Optimizing Binary Decision Diagrams for Interpretable Machine Learning Classification[END_REF] is the only exact method of learning optimal BDDs in the context of ML. The authors proposed a SAT model to learn optimal BDDs with the smallest sizes that correctly classify all examples. In their approach, the depth of the BDD is not restrained. In fact, it is possible that the constructed BDD is small in size (number of nodes) and high in depth. As the BDD is ordered, this approach could not limit the number of features used, making it not quite comparable with our proposition. Another related work is in [START_REF] Hu | Learning Optimal Decision Trees with MaxSAT and its Integration in AdaBoost[END_REF] where the authors consider a MaxSAT model to learn optimal decision trees minimizing the classification error within a limited depth. The usage of the same solving methodology with the same objective function and the depth limit, makes these two MaxSAT models comparable. Finally, in order to increase the scalability of our approach, we propose a heuristic extension based on a simple pre-processing step. For the sake of space, some details are left in a technical report [START_REF] Hu | Optimizing Binary Decision Diagrams with MaxSAT for classification[END_REF].

Technical Background Classification

Consider a dataset E = {e q , . . . , e M } with M examples.

Each example e q ∈ E is characterized by a list of binary features L q = [f 1 , . . . , f K ] and a binary target cl q , representing the class of the example (cl q ∈ {0, 1}). The data set is partitioned into E + and E -, where E + (respectfully E -) is the set of positive (respectifully negative) examples. That is, cl q = 1 iff e q ∈ E + and cl q = 0 iff e q ∈ E -. We assume that, ∀1 ≤ q, q ′ ≤ M , L q = L q ′ implies cl q = cl q ′ . Let ϕ be the function defined by ϕ(L q ) = cl q , ∀q ∈ [1, M ]. The classification problem is to compute a function γ (called a classifier) that matches as accurately as possible the function ϕ on examples e q of the training data and generalizes well on unseen test data.

Binary Decision Diagrams

Binary Decision Diagrams (BDDs) are used to provide compact representation of Boolean functions. Let [x 1 , . . . , x n ] be a sequence of of n Boolean variables. A BDD is a rooted, directed, acyclic graph G. The vertex set V of G contains two types of vertices. A terminal vertex v is associated to a binary value: value(v) ∈ {0, 1}. A nonterminal vertex v, is associated to a Boolean variable x i and has two children left(v), right(v) ∈ V. In this case, index(v) = i ∈ {1, . . . , n} is the index of the Boolean variable associated to v.

We assume that all BDDs are ordered and reduced. These two restrictions are widely considered in the literature as they guarantee a unique BDD for a given Boolean function. The restriction ordered indicates that for any nonterminal vertex v, index(v) < index(left(v)) and index(v) < index(right(v)). The restriction reduced indicates that the graph contains no nonterminal vertex v with left(v) = right(v), nor does it contain distinct nonterminal vertices v and v ′ having isomorphic rooted sub-graphs. Therefore, given an ordered reduced BDD G with root v, the associated Boolean function can be recursively obtained with the Shannon expansion process [START_REF] Shannon | A symbolic analysis of relay and switching circuits[END_REF].

Let g be a Boolean function defined over a sequence X = [x 1 , . . . , x n ] of n Boolean variables. The function g can be represented by a truth table that lists the 2 n values of all assignments of the n variables. The value of the truth table is therefore associated to a string of 2 n binary values. A truth table β of length 2 n is said to be of order n. A truth table β of order n > 0 has the form β 0 β 1 , where β 0 and β 1 are truth tables of order n -1, and β 0 and β 1 are called subtables of β. The subtables of subtables are also considered to be subtables, and a table is considered as a subtable of itself. A bead of order n is a truth table T of order n that does not have the form αα where α is a subtable of T . The beads of g are the subtables of its truth table that happen to be beads. Proposition 1 from [START_REF] Knuth | The Art of Computer Programming, Volume 4, Fascicle 1: Bitwise Tricks & Techniques; Binary Decision Diagrams[END_REF] relates truth table and binary decision diagram for the same Boolean function. Proposition 1. All vertices in V of a binary decision diagram G, are in one-to-one correspondence with the beads of the Boolean function g it represents.

Based on Proposition 1, we can produce the ordered and reduced binary decision diagram of a Boolean function by finding its beads and combine its beads with its sequence of variables. Example 1. Consider the Boolean function from [START_REF] Knuth | The Art of Computer Programming, Volume 4, Fascicle 1: Bitwise Tricks & Techniques; Binary Decision Diagrams[END_REF]:

g 1 (x 1 , x 2 , x 3 ) = (x 1 ∨ x 2 ) ∧ (x 2 ∨ x 3 ) ∧ (x 1 ∨ x 3 ).
The binary string associated to its truth table β is 00010111. The beads of β are {00010111, 0001, 0111, 01, 0, 1}.

From Proposition 1, we can draw the BDD with the beads found, shown as the left part of Figure 1. The dashed (solid) line of each vertex indicates the left (right) child. Then, we can replace the beads by vertices associated with the sequence of Boolean variables. The final binary decision diagram for g 1 is shown as the right part of Figure 1.

Oblivious Read-Once Decision Graphs

Oblivious Read-Once Decision Graphs (OODGs) are proposed in [START_REF] Kohavi | Bottom-Up Induction of Oblivious Read-Once Decision Graphs[END_REF] to overcome some limitations of decision trees for multi-classification, like replication and fragmentation problem. We refer the readers to [START_REF] Kohavi | Oblivious Decision Trees, Graphs, and Top-Down Pruning[END_REF][START_REF] Kohavi | Bottom-Up Induction of Oblivious Read-Once Decision Graphs[END_REF] for details on OODGs. An OODG is a rooted, directed, acyclic graph, which contains terminal category nodes labelled with classes to make decisions, and non-terminal branching nodes labelled with features to make splits. The property "read-once" indicates that each feature occurs at most once along any path from the root to a category node. The property "levelled" indicates that the nodes are partitioned into a sequence of pairwise disjoint sets, representing the levels, such that outgoing edges from each level terminate at the next level. The property "oblivious" extends the idea of "levelled" by guaranteeing that all nodes at a given level are labelled by the same feature.

For the classification process, top-down and bottom-up heuristic methods for building OODGs are proposed in [START_REF] Kohavi | Oblivious Decision Trees, Graphs, and Top-Down Pruning[END_REF][START_REF] Kohavi | Bottom-Up Induction of Oblivious Read-Once Decision Graphs[END_REF]. Here, we introduce briefly the top-down heuristic method, which is similar to the heuristic methods C4.5 and CART for computing decision trees. The top-down heuristic induction for OODG with given depth contains three critical phases: (1) selecting a sequence of features with the help of mutual information (the difference of conditional entropy (Cover and Thomas 2006)); ( 2) growing an oblivious decision tree (ODT) by splitting the dataset with features in the sequence selected; and (3) merging isomorphic and compatible subtrees from top to down to build the OODG. When building the ODT, the algorithm marks nodes that capture no example of the dataset as "unknwon". For the merging phase, two subtrees are compatible if at least one root is labelled as "unknown", or if the two root nodes are labelled with same feature and their corresponding children are the roots of compatible subtrees. The ODT grown could make classifications directly by assigning "unknown" nodes with the majority class of their parents. Figure 2 shows an example of two compatible subtrees and the merged tree, where "unknown" nodes are labelled as "u". Merging compatible subtrees changes the bias by assuming that a "unknown" node is likely to behave the same as another child if they belong to compatible subtrees.

f 1 u f 2 1 u 1 y e s n o v 1 v 2 v 3 f 1 f 2 f 2 1 1 1 1 0 u y e s n o v 1 v 2 v 3 v 1 v 2 v3 f 1 f 2 f 2 1 1 1 1 0 1 y e s n o v 1 v 2 v 3 v 1 v 2 v 3
In binary classification for binary datasets, OODGs could be considered equivalent to BDDs, as the properties "oblivious" and "read-once" for OODGs are same as property "ordered" for BDDs. In addition, the use of merging compatible subtrees could also be applied for BDDs.

SAT and MaxSAT

We use standard terminology for Boolean Satisfiabily (Biere et al. 2009). A literal is a Boolean variable or its negation, and a clause is a disjunction of literals. An assignment of variables satisfies a clause if one of its literals is true. Given a set of Boolean variables and a set of clauses defined over these variables, the SAT problem can be defined as finding an assignment of the variables such that all the clauses are satisfied. Maximum Satisfiability (MaxSAT) is an optimization version of the SAT problem, where the clauses are partitioned into hard and soft clauses. Here we consider the Partial MaxSAT problem, that is to find an assignment of the Boolean variables that satisfies all the hard clauses and maximizes the number of satisfied soft clauses.

(Max)SAT-Based Model for Binary Decision Diagrams

In this section, we present our approach for learning BDDs for binary classification using SAT and MaxSAT.

Problem Definition

We firstly consider the following decision problem for classification with BDD in a given depth. To realize the SAT encoding, we introduce two sets of Boolean variables as follow:

• a i r : the variable a i r is 1 iff feature f r is selected as i-th feature in the feature ordering, where i = 1, . . . , H, r = 1, . . . , K.

• c j : the variable c j is 1 iff the j-th value of the truth table is 1, where j = 1, . . . , 2 H . The set of variables a i r guarantees the ordered restriction. Then, we introduce two constraints (1) and ( 2) for the feature ordering. Constraint 1 ensures that any feature f r can be selected at most once.

H i=1 a i r ≤ 1, r = 1, . . . , K (1)
Then, there is exacty one feature selected for each index of the feature ordering.

K r=1 a i r = 1, i = 1, . . . , H (2)
We use the classical sequential counter encoding proposed in [START_REF] Sinz | Towards an Optimal CNF Encoding of Boolean Cardinality Constraints[END_REF] to model constraints (1) and (2) as a Boolean formula.

The truth table we are looking for is the binary string of the values of variables c 1 c 2 . . . c 2 H . To avoid the first feature selected makes useless split, we need to make sure that the truth table is a bead.

2 H-1 j=1 (c j ⊕ c j+2 H-1 ) (3)
There is a relationship between the values of a truth table and the assignments of the given sequence of Boolean variables. For example, the first value of a truth table corresponds to the assignment that x 1 = 0 and x 2 = 0. Therefore, we define the following function to obtain the value of the i-th feature in the feature ordering of size H given the j-th value in the truth table.

rel(i, j) = ⌊ j -1 2 H-i ⌋ mod 2, i ∈ [1, H], j ∈ [1, 2 H ] (
4) For an example e q ∈ E, we denote the value of the feature f r as σ(r, q). If rel(i, j) = σ(r, q), it indicates that for example e q , the feature f r can be at the i-th position in the feature ordering to produce the j-th value in the truth table.

To classify all examples correctly, we ensure that no example follows an assignment in the truth table leading to its opposite class. Thus, we propose the following constraints for classification. Let e q ∈ E + , for all j = 1, . . . , 2 H :

¬c j → H i=1 K r=1 (a i r ∧ rel(i, j) ⊕ σ(r, q)) (5)
That is, for every positive example e q , any variable c j assigned to 0 must be associated to an assignment of features that contains at least one feature-value that is not coherent with e q . For negative examples, we use a similar idea. Let e q ∈ E -, for all j = 1, . . . , 2 H : The two sets of variables are: {a 1 1 , a 2 1 , a 1 2 , a 2 2 , a 1 3 , a 2 3 , a 1 4 , a 2 4 }, and {c 1 , c 2 , c 3 , c 4 }. The constraints 1, 2, and 3 are:

c j → H i=1 K r=1 (a i r ∧ rel(i, j) ⊕ σ(r, q)) (6) E0 f1 f2 f3 f4 c e1 1 0 1 0 0 e2 1 0 0 1 0 e3 0 0 1 0 1 e4 1 1 0 0 0 e5 0 0 0 1 1 e6 1 1 1 1 0 e7 0 1 1 0 0 e8 0 0 1 1 1 Table 1: A binary classification dataset f 1 f 2 0 1 0
a 1 1 + a 2 1 ≤ 1, a 1 2 + a 2 2 ≤ 1, a 1 3 + a 2 3 ≤ 1, a 1 4 + a 2 4 ≤ 1 a 1 1 + a 1 2 + a 1 3 + a 1 4 = 1, a 2 1 + a 2 2 + a 2 3 + a 2 4 = 1 (c 1 ⊕c 3 ) ∨ (c 2 ⊕ c 4 )
For classification constraints (i.e., 5 and 6), we show the encoding of e 1 ∈ E -with for value c 1 . The encoding for other examples and other values is similar.

c 1 → (a 1 1 ∧ 0 ⊕ 1) ∨ (a 1 2 ∧ 0 ⊕ 0) ∨ (a 1 3 ∧ 0 ⊕ 1) ∨(a 1 4 ∧ 0 ⊕ 0) ∨ (a 2 1 ∧ 0 ⊕ 1) ∨ (a 2 2 ∧ 0 ⊕ 0) ∨(a 2 3 ∧ 0 ⊕ 1) ∨ (a 2 4 ∧ 0 ⊕ 0)
This could be simplified as follow: 2 illustrates the relationship between the values of truth table and the assignments of the given sequence Boolean variable of size 2. Figure 4 shows the corresponding BDD. This BDD classifies all examples of the dataset E 0 correctly, also provides more compact representation than the decision tree shown in Figure 3.

¬c 1 ∨ a 1 1 ∨ a 1 3 ∨ a 2 1 ∨ a 2 3 x1 = f1 x 2 = f 2 0 0 c 1 = 1 0 1 c 2 = 0 1 0 c 3 = 0 1 1 c 4 = 0
We refer to this first SAT encoding for P bdd (E, H) as BDD1. The size of BDD1 is given in Proposition 2.

Proposition 2. For a P bdd (E, H) problem with K binary features and M examples, the encoding size (in terms of the number of literals used in the different clauses

) of BDD1 is O(M × H × K × 2 H ).
Proof. Notice first that j ranges from 1 to 2 H , i ranges from 1 to H, and r ranges from 1 to K. The term M × 2 H results from constraint (5) and ( 6), each contains O(H × K) literals. For the remaining constraints, it is O(H × K) for constraints (1) and (2), O(2 H ) for constraint [START_REF] Breiman | Classification and Regression Trees[END_REF].

The size of BDD1 is quite huge due to the size of clauses generated by constraints ( 5) and ( 6) for classification. This makes BDD1 impractical in practice.

An Improved SAT Model for P bdd (E, H)

In order to reduce the size of BDD1, we propose new classification constraints to replace constraints ( 5) and ( 6). The idea is that every positive (respectivery negative) example follows an assignment leading to a positive (respectively negative) value of the truth table. We introduce a new set of Boolean variables:

• d q i : The variable d q i is 1 iff for example e q the value of the i-th feature selected in feature ordering is 1, where i = 1, . . . , H, q = 1, . . . , M .

Then, We describe constraints that relate the values of features for each example e q ∈ E, for i = 1, . . . , H, r = 1, . . . , K:

a i r → d q i if σ(q, r) = 1 a i r → ¬d q i if σ(q, r) = 0 (7)
Let e q ∈ E + , we have 2 H constraints for classifying examples correctly:

¬d q 1 ∧ ¬d q 2 ∧ • • • ∧ ¬d q H-1 ∧ ¬d q H → c 1 ¬d q 1 ∧ ¬d q 2 ∧ • • • ∧ ¬d q H-1 ∧ d q H → c 2 . . . d q 1 ∧ d q 2 ∧ • • • ∧ d q H-1 ∧ d q H → c 2 H (8)
That is, any positive example follows an assignment of the feature ordering that leads to a positive value in the truth table.

Similarly, for any e q ∈ E -, we also have 2 H constraints:

¬d q 1 ∧ ¬d q 2 ∧ • • • ∧ ¬d q H-1 ∧ ¬d q H → ¬c 1 ¬d q 1 ∧ ¬d q 2 ∧ • • • ∧ ¬d q H-1 ∧ d q H → ¬c 2 . . . d q 1 ∧ d q 2 ∧ • • • ∧ d q H-1 ∧ d q H → ¬c 2 H (9)
We refer to this new SAT encoding for P bdd (E, H) as BDD2. The encoding size of BDD2 is given in Proposition 3. Proposition 3. For a P dd (E, H) problem with K binary features and M examples, the encoding size of the SAT encoding

(BDD2) is O(M × H × (2 H + K)).
Proof. The term M × H × K results from constraint (7). For constraints ( 8) and (9), for each example, there are 2 H clauses containing H + 1 literals. The term M × H × 2 H results from that.

Propositions 2 and 3 show a clear theoretical advantage of BDD2 compared to BDD1 in terms of the encoding size, thus scalability.

MaxSAT Model for P *

bdd (E, H) :

We now present a MaxSAT encoding for the optimization problem P * bdd (E, H). That is, given a set of examples E, find a binary decision diagram of depth H that maximises the number of examples correctly classified.

We transform the SAT encoding of BDDs into a MaxSAT encoding following a simple technique. The idea is to keep structural constraints as hard clauses and classification constraints as soft clauses. We consider BDD2 as it has a reduced size. Constraints (1), ( 2), ( 3) and ( 7) are kept as hard clauses. To classify the examples, we declare all clauses of constraints ( 8) and ( 9) as soft clauses. For any example e q , the number of satisfied soft clauses associated to e q is either 2 H (indicating e q is classified correctly), or 2 H -1 (indicating e q is classified wrongly). Therefore, the objective of maximising the number of satisfied soft clauses is equivalent to maximise the number of examples that are correctly classified.

Merging Compatible Subtrees

Consider a BDD G found by a MaxSAT solver and its associated truth table β. Based on the feature ordering of G, it is possible that some values in β capture no (training) example (Equivalent to "unknown" nodes for OODG). Such values are decided by the MaxSAT solver in an arbitrary way, which gives a certain bias in generalisation. We propose to merge compatible subtrees in G in order to handle this bias. This will result in changing some values in the truth table β (i.e. the arbitrary ones decided by MaxSAT).

We propose a post-processing procedure to merge compatible subtrees using the following three phase: (1) update the truth table β by replacing the values of β that capture no examples with a special value "u"; (2) for each level, check the beads, where "u" can be used to match 1 or 0, and create a node for each bead; (3) for each level, after creating the nodes, check the matches between all subtables of the next level. For matched subtables, update the corresponding beads of current level to eliminate the "u" values.

Experimental Results

We present our large experimental study to evaluate empirically our propositions on different aspects2 . We consider datasets from CP4IM 3 . These datasets are binarized with the one-hot encoding. Preliminary experiments on the decision problem P bdd (E, H) confirm the great improvements the encoding size of BDD2 compared to BDD1, as shown in propositions 2 and 3. However, solving the optimization problem of finding a BDD that classifies all examples in the dataset correctly with a minimum depth using linear calls to BDD1 and BDD2 was hard. We therefore focus on the optimisation problem P * bdd (E, H). This is also motivated by the fact that classifying all examples correctly induces overfitting.

At first, we evaluate the prediction performance between the proposed MaxSAT-BDD model and the heuristic method, ODT and OODG. Next, we compare our model with an exact method for building decision trees using MaxSAT [START_REF] Hu | Learning Optimal Decision Trees with MaxSAT and its Integration in AdaBoost[END_REF]) in terms of prediction quality, model size, and encoding size. Finally, we propose and evaluate a simple heuristic version of our encoding to tackle scalability. For each dataset, we use random 5-fold cross-validation with 5 different seeds. All experiments were run on a cluster using Xeon E5-2695 v3@2.30GHz CPU and running xUbuntu 16.04.6 LTS. The MaxSAT solver we used is Loandra [START_REF] Berg | Core-Boosted Linear Search for Incomplete MaxSAT[END_REF], an efficient incomplete MaxSAT solver that return the best solution found within a limited computation time or report optimality. For each experiment, the time limit for generating formulas and the time limit for solver are set to 15 minutes.

Comparison with Existing Heuristic Approaches

We consider the P * bdd (E, H) problem with 5 different depths H ∈ {2, 3, 4, 5, 6}. We compare our MaxSAT-BDD model with the approach proposed in [START_REF] Kohavi | Oblivious Decision Trees, Graphs, and Top-Down Pruning[END_REF] to learn ODT and OODG. For the heuristic methods, as described in the background section, after merging isomorphic and compatible subtrees of ODT, the corresponding OODG changes the bias for those "unknown" nodes. In fact, different bias affects the prediction for unseen examples, but not for learned examples. Therefore, the training accuracy of ODT and OODG are equals, but the testing accuracy could be different. This fact also suits to the MaxSAT-BDD with different biases for the "unknown" nodes. In this experiment, we consider three different biases: assigning for each unknown node the majority class of its branch (denoted as MaxSAT BDD-P), merging compatible subtrees (MaxSAT BDD-C), and the class decided by the MaxSAT solver (MaxSAT BDD-S).

The left scatter plot in Figure 5 presents the comparison of the average training accuracy between OODG and MaxSAT-BDD model. In this figure, different datasets are marked with different colors, and different depths are labelled with points of different sizes. From the scatter plot, we observe that the average training accuracy of both approaches increase with the increase of depth. Overall, the MaxSAT-BDD model performs better than the heuristic OODG in training accuracy.

The middle boxplots of Figure 5 show the average testing accuracy of MaxSAT-BDD with different biases, ODT, and OODG using different depths averaged over all datasets. The white line and green triangle of each box indicate the median and the average value, respectively. Several observations from the boxplots are presented as follows. At first, for each bias, increasing the depth could improve the prediction performance. However, compared to ODT and OODG, all biases chosen for MaxSAT-BDD get less improvements with increasing depths. Next, there are slight differences for different biases of MaxSAT-BDD, indicating that it is quite robust. Then, generally, MaxSAT-BDD gets better prediction performance than ODT and OODG, in particular when depths are small. We noticed also that when the depth is 2, all datasets (except one), MaxSAT-BDD reports optimality.

Comparison with an Exact Decision Tree Approach

The purpose of this experiment is to compare our proposition with the exact method for learning decision trees using the same solving approach (MaxSAT). For MaxSAT-BDD, we consider only the bias of merging compatible subtrees (MaxSAT BDD-C) since no substantial difference was observed. We consider different values for depth: H ∈ {2, 3, 4, 5, 6}. For MaxSAT-BDD, the depth also corresponds to the number of selected features, whereas for MaxSAT-DT the depth chosen here indicate maximum depth of the BDD.

Table 3 presents the results of evaluation. In the column "Dataset", the dataset size (left) and the number of binary feature (right) are shown under the name. The column "Size" and "E S" indicate the number of nodes of model and the encoding size (number of literals in 100 thousands). The best values are marked in bold.

The results in Table 3 show that the MaxSAT-BDD is competitive to MaxSAT-DT in terms of prediction quality. In most cases, the training and testing accuracy of these two approaches are close. However, the size of the models are always smaller with MaxSAT-BDD. The difference grows bigger when the depth increases. The reduction in model size provides better intrepretability. Moreover, sometimes, compared to the optimal BDDs found via MaxSAT-BDD, the optimal decision trees found via MaxSAT-DT make useless splits improving no prediction performance, as we can see in the case of the datasets "car" and "hypothyroid" with depth 2. We observe also that MaxSAT-BDD has always a (far) lighter encoding size than MaxSAT-DT which makes it easier to handle and to report optimality in limited time.

Evaluation of a Heuristic MaxSAT-BDD Method

To increase the scalability of our model, we propose a heuristic version of MaxSAT-BDD using CART [START_REF] Breiman | Classification and Regression Trees[END_REF], an efficient and scalable heuristic for learning decision trees. The idea is to run CART as a pre-processing step in order to choose a subset of (important) features. Then apply the MaxSAT-BDD approach using only the selected subset of features.

We first evaluate the prediction quality of this heuristic method with the original MaxSAT-BDD. The size of the encoding is in favor of the heuristic approach as expected, making large problems treatable. The results of average testing accuracy are shown in the right scatter of Figure 5. This heuristic approach is clearly very competitive to the original MaxSAT-BDD in generalization performance. Especially, for those datasets with more features, the heuristic approach obtains better prediction performance than the original within the same limited time. As a consequence, compared to OODG, our heuristic approach performs much better in terms of prediction quality. We also observed no significant prediction quality difference with CART.

Conclusion

We propose exact and heuristic methods for optimizing binary decision diagrams (BDDs) based on the (Maximum) Boolean Satisfiability framework. Our large experimental studies show clear benefits of the proposed framework in terms of prediction quality and interpretability compared to the existing heuristic approach. Besides, our approach has competitive prediction performance with far lighter encoding size compared to a similar approach for building decition trees.

In the future, it would be interesting to extend the proposed approach for multi-valued classification. Moreover, a deeper investigation of BDDs with other interpretable models (such as decision rules and decision sets) will be beneficial since it enhances alternatives for explainable AI.

Figure 1 :

 1 Figure 1: The Binary decision diagram for g 1 (x 1 , x 2 , x 3 )

Figure 2 :

 2 Figure 2: An example of two compatible subtrees (the left two) and the merged tree (the right one) from[START_REF] Kohavi | Oblivious Decision Trees, Graphs, and Top-Down Pruning[END_REF]

Figure 3 :

 3 Figure 3: Decision Tree found

Figure 5 :

 5 Figure 5: The left scatter shows the average training accuracy of OODG and MaxSAT model. The middle boxplots show the average testing accuracy with different biases: MaxSAT BDD-P(blue), MaxSAT BDD-C(yellow), MaxSAT BDD-S(green), ODT(red), OODG(purple). The right scatter shows the average testing accuracy of MaxSAT-BDD and its heuristic approach.

  Notice that the algorithm for P bdd (E, H) can be used to the alternative problem of optimizing a BDD that classifies all examples in the dataset correctly with a minimum depth. For that purpose, one can use a linear search that takes an initial depth H 0 as input and progressively increases or decreases this value depending on the result of solving P bdd (E, H). Given a set of examples E, find a BDD of depth H that maximises the number of examples in E that are correctly classified.We propose an initial SAT model for the decision problem P bdd (E, H). Then, we propose an improved version in tighter formula size. Finally, we show how the improved SAT model for P bdd (E, H) can be used effectively to solve the optimization problem P * bdd (E, H) with MaxSAT.

	Next, we consider another optimization problem for the
	classification with BDD in a limited depth.
	• P * bdd (E, H) :

• P bdd (E, H) : Given a set of examples E, is there a BDD of depth H that classifies correctly all examples in E? SAT Model for P bdd (E, H) As shown before, a BDD of depth H could be generated from the combination of a sequence of Boolean variables of size H: [x 1 , . . . , x H ], and a truth table of order H associated to a Boolean function. To solve the classification problem P bdd (E, H), we then have to find a sequence of binary features of size H that maps one-to-one the sequence of Boolean variables, and a truth table associated to a Boolean function that well-classified all examples. We denote the sequence of binary features found as feature ordering. Therefore, the SAT encoding consists of two parts: • Part 1: Constraints for selecting features of the dataset into the feature ordering of size H. • Part 2: Constraints for generating a truth table that classifies all examples of E correctly with the selected feature ordering.

Table 2 :

 2 Truth table solution forThe values of truth table found by the SAT model are shown in Table2, the feature ordering is [f 1 , f 2 ]. Moreover, Table

			f 1
		f 2
		1	0
	BDD of depth 2 classifying all ex-ample of E 0	Figure 4: The BDD found

  .49 8.64 0.9 87.89 82.85 20.92 7.64 5 87.51 85.83 13.32 1.38 90.1 79.34 47.80 23.91 6 88.57 81.12 19.84 2.24 90.34 81.31 37.32 83.86 soybean (630/50) 2 90.48 90.48 4 1.08 91.27 91.27 7 2.56 3 91.39 90.41 6.52 1.7 95.45 94.7 15 6.23 4 93.24 93.21 9.04 2.45 97.25 95.9 22.20 16.02 5 94.31 92.95 11.92 3.53 97.96 95.3 40.60 45.53 6 96.07 95.52 14.88 5.34 98.27 96.03 33.40 145.99 splice-1 (3190/287) 2 84.04 84.04 4 29.66 84.22 83.17 6.92 55.52 3 87.25 86.94 5.44 44.9 87.79 87.37 11.32 123.16 4 88.3 88.04 7.24 60.83 86.52 85.64 16.60 271.79 5 71.99 70.53 10.28 78.39 77.37 76.32 21.88 622.68 6 62.92 61.89 16.28 99.63 60.36 58.95 29.40 1540.61 vote (435/48) 2 95.68 95.22 3.76 0.72 96.21 95.03 7 1.83 3 96.69 94.57 5.56 1.14 97.39 93.79 13.96 4.66 4 97.40 94.39 8.16 1.65 98.62 94.57 21.16 12.65 5 98.21 94.57 12.4 2.38 99.47 93.84 30.52 38.2 6 98.93 93.98 18.44 3.62 99.62 94.76 35.40 129.24

	Datasets H	MaxSAT BDD-C Train Test Size E S Train Test Size E S MaxSAT-DT
	anneal (812/89) audiology (216/146) australian (653/124) cancer (683/89)	2 82.92 82.19 5 3 84 83.55 7 4 84.58 83.84 9.4 5.21 86.05 84.78 18.68 31.55 2.41 83.18 82.14 6.84 5.27 3.72 85.07 84.66 12.68 12.62 5 85.33 83.92 11.72 7.11 86.44 84.88 23.88 86.53 6 86.26 83.70 14.68 9.95 87.6 85.76 39.16 266.67 2 94.91 94.92 4 1.06 95.49 94.92 7 3.14 3 96.78 95.84 5.04 1.64 97.82 95.56 11.56 8.88 4 97.73 95.56 6.96 2.26 99.51 94.54 19.08 27.21 5 98.40 94.44 9.88 2.98 99.95 93.98 27 91.53 6 99.17 95.84 14.28 3.96 99.86 94.08 24.12 332.36 2 86.70 85.94 4.72 2.68 86.93 85.33 6.68 5.97 3 87.45 84.81 5.32 4.11 88.09 84.87 13.08 14.61 4 88.45 86.03 7.4 5.68 88.74 85.18 17.48 37.76 5 89.36 85.91 10.44 7.59 89.28 84.75 22.52 107.64 6 90.05 85.7 17.32 10.25 89.49 84.84 27.08 343.36 2 93.88 93.59 4 2.03 94.91 94.2 7 4.56 3 95.02 93.91 5.84 3.14 96.6 94.73 15 11.09 4 96.06 95.49 7.96 4.39 97.34 94.17 21 28.38 5 95.94 93.74 10.68 5.99 97.99 94.35 29.32 80.09 6 96.84 94.35 14.8 8.38 98.87 93.41 45.72 253.69 2 85.53 85.53 4 1.33 85.53 85.53 6.84 3.2	MaxSAT BDD-C Train Test Size E S Train Test Size E S MaxSAT-DT 2 79.04 72.57 4 0.95 80.76 72.84 7 Datasets H 2.56 3 85.07 83.37 6 1.47 85.68 76.55 12.84 6.89 cleveland 4 86.32 79.46 7.84 2.06 86.77 76.75 17.80 20.08 (296/95) 5 88.65 78.72 13.08 2.79 87.26 74.45 23.96 64.68 6 90.74 77.29 21.04 3.87 88.58 75.81 28.84 228.48 2 77.83 77.01 4 7.79 86.92 86.92 7 15.51 3 90.43 90.43 5.28 12.05 93.81 93.79 12.44 34.3 kr-vs-kp 4 94.09 94.09 7.56 17.03 94.32 94.14 17.24 75.38 (3196/73) 5 94.34 94.18 9.52 23.64 94.85 94.69 25.40 171.71 6 92.80 92.55 11.52 33.94 93.91 93.69 29.32 422.77 2 84.46 83.23 4 0.35 86.01 79.27 7 1.23 3 86.76 78.35 5.92 0.56 91.93 80.54 14.68 3.67 lymph 4 90.54 82.4 8.72 0.79 94.56 78.46 20.20 11.79 (148/68) 5 93.51 83.6 13.52 1.09 97.09 82.46 27.08 41.31 6 95.88 84.82 17.64 1.57 99.59 80.92 46.60 155.03 2 82.80 81.6 4 0.37 82.92 81.01 6.76 1.05 3 83.84 80.43 5.3 0.6 86.16 82.97 13.88 2.72 tumor (336/31) 4 85.52 82
	car (1728/21)	3 88.40 87.41 5.08 2.2 89.25 87.45 12.68 7.18 4 89.84 88.54 6.84 3.44 91.62 89.68 20.36 16.25 5 91.13 89.91 9.6 5.58 93.78 92.77 29.56 39.0
		6 93.51 92.99 13.36 9.71 95.8 95.06 31.96 104.45
		2 97.84 97.84 4	9.27 97.84 97.84 5.96 18.22
	hypothyroid (3247/86)	3 98.09 98.04 5.12 14.28 98.14 97.82 9.72 40.3 4 98.27 98.13 6.72 20.01 98.38 98.01 15.40 88.55 5 98.30 98.05 9.28 27.4 98.45 98 20.04 201.63
		6 98.37 97.95 13.68 38.54 98.46 97.91 33.16 495.76
		2 95.13 95.13 4 29.92 96.9 96.9 7	56.53
	mushroom (8124/112)	3 97.74 97.77 6.8 45.81 99.9 99.9 13.72 122.72 4 98.78 98.74 9 63.51 100 100 19.80 260.39 5 98.63 98.57 11.32 85.37 100 100 23.40 557.11
		6 97.28 97.10 14.6 116.59 100 100 27.56 1237.69
		2 71.05 68.35 4	0.93 71.1 67.49 5.96 2.23
	tic-tac-toe (958/27)	3 74.91 72.36 6.16 1.5 77.15 73.55 11.48 5.2 4 76.87 74.22 8.84 2.29 82.47 78.68 20.60 12.51 5 81.86 80.31 13.88 3.57 83.08 79.50 28.44 32.83
		6 84.82 80.08 24.16 5.95 84.25 80.86 38.12 97.95

Table 3 :

 3 Comparison Results between MaxSAT-BDD and MaxSAT-DT.

The source code and the datasets are available online at https://gitlab.laas.fr/hhu/bddencoding

https://dtai.cs.kuleuven.be/CP4IM/datasets/


Karol Tylmann
email: k.tylmann@ug.edu.pl

Vincent R Rinterknecht

Piotr P Woźniak

Didier Bourles

Irene Schimmelpfennig

Valery Guillou

Georges Aumaitre

Karim Keddadouche

Didier Bourlès

Valery Guilou

Aster Team

Georges Aumaître

The Local Last Glacial Maximum of the southern Scandinavian Ice Sheet front: Cosmogenic nuclide dating of erratics in northern Poland

Keywords: Local Last Glacial Maximum, Leszno (Brandenburg) Phase, Poznań (Frankfurt Phase), cosmogenic nuclide dating, Scandinavian Ice Sheet, erratics

Introduction

The Last Glacial Maximum (LGM) is defined as the most recent period with the maximum global integrated volume of ice sheets and the corresponding minimum global sea level [START_REF] Mix | Environmental processes of the ice age: land, oceans, glaciers (EPILOG)[END_REF][START_REF] Clark | The Last Glacial Maximum[END_REF][START_REF] Hughes | Timing of glaciation during the last glacial cycle: evaluating the concept of a global 'Last Glacial Maximum' (LGM)[END_REF] and dated to 23-19 cal ka BP or 24-18 cal ka BP (cf. [START_REF] Mix | Environmental processes of the ice age: land, oceans, glaciers (EPILOG)[END_REF][START_REF] Hughes | A stratigraphical basis for the Last Glacial Maximum (LGM)[END_REF]. The termination of the

LGM climatic cooling resulted in the progressive shrinking of the inland ice, which for most of the Northern Hemisphere ice sheets is dated at ~20-17 ka [START_REF] Hughes | A stratigraphical basis for the Last Glacial Maximum (LGM)[END_REF]. However, numerical ages (based on calibrated 14 C, terrestrial cosmogenic nuclides -TCN or luminescence dating) of terrestrial deposits show significant regional variability in the timing of the ice sheets' maximum extent and the commencement of their retreat. Most of the ice sheets reached their maximum extents in a 7.0 ka time interval between 26.5 ka and 19.0 ka [START_REF] Clark | The Last Glacial Maximum[END_REF], but the maximum extents of specific ice sheets (Local Last Glacial Maxima -LLGM) were found to be asynchronous and not always in agreement with the global LGM [START_REF] Hughes | Timing of glaciation during the last glacial cycle: evaluating the concept of a global 'Last Glacial Maximum' (LGM)[END_REF].

The Scandinavian Ice Sheet (SIS) was one of the main Northern Hemisphere ice sheet and covered a significant part of northern Europe during the LGM. Timing of the last SIS maximum expansion and retreat is important to understand climatic fluctuations at the end of the Marine Isotope Stage (MIS) 2. The LLGM to the south of the Baltic Basin (Fig. 1A) is dated at ~24 ka in Germany and western Poland, and at ~19 ka in eastern Poland based on calibrated 14 C and OSL ages [START_REF] Marks | Timing of the Late Vistulian (Weichselian) glacial phases in Poland[END_REF] after [START_REF] Stankowska | Maximum extent of the Vistulian ice sheet in the vicinity of Konin, Poland: a geomorphological, sedimentological and radiometric evidence[END_REF][START_REF] Wysota | Record of the Vistula ice lobe advances in the Late Weichselian glacial sequence in north-central Poland[END_REF] and at 21-19 ka in Belarus and Lithuania based on calibrated 14 C and 10 Be ages [START_REF] Rinterknecht | Timing of the Last Deglaciation in Belarus[END_REF][START_REF] Rinterknecht | Timing of the last deglaciation in Lithuania[END_REF]. Thus, the maximum extent of the of the last SIS in Poland was probably asynchronous with earlier maximum expansion of the ice sheet in western and central Poland and subsequent maximum expansion further to the east (Fig. 1A andB).

Nevertheless, there is still a relatively small number of geochronological data constraining the timing of LLGM in north Poland. Specifically, there is a lack of direct dating (TCN) of moraines correlated with the maximum limit of the last ice sheet in this region.

Here we present a new set of TCN exposure ages of erratic boulders located within the southern front of the last SIS in northern Poland. It provides the first exposure ages of erratics located on the pre-LLGM moraines and new exposure ages on the LLGM moraines in Poland.

Additionally, published 14 C ages were calibrated with the most recent radiocarbon calibration curve and used as a background indicator of the possible time window for the LLGM. We discuss our results in a wider chronological context including specifically the growing number of luminescence ages available for the southern fringe of the last SIS (e.g. Houmark-Nielsen and Kjaer, 2003; [START_REF] Wysota | Record of the Vistula ice lobe advances in the Late Weichselian glacial sequence in north-central Poland[END_REF][START_REF] Lüthgens | Age of the Pomeranian ice-marginal position in northeastern Germany determined by Optically Stimulated Luminescence (OSL) dating of glaciofluvial sediments[END_REF][START_REF] Bitinas | New sights into the last deglaciation of the southeastern flank of the Scandinavian Ice Sheet[END_REF][START_REF] Lasberg | Chronology of Late Weichselian glaciation in the western part of the East European Plain[END_REF], as well as existing reconstructions of the SIS evolution (e.g. [START_REF] Hughes | The last Eurasian ice sheets e a chronological database and time-slice reconstruction[END_REF][START_REF] Stroeven | Deglaciation of Fennoscandia[END_REF]. The objective of this paper is to date the asynchronous maximum extent of the last SIS in northern Poland based on new exposure ages, and on the extensive dataset of TCN exposure ages already published.

Study area

The study area is located in NW Poland and covers ~20% (~500 km) of the LLGM ice sheet margin south of the Baltic Sea. The investigated region is situated within the eastern part of the extensive, lobe-shaped fringe of the last SIS (Fig. 1A). The region is large enough to find a considerable number of erratics suitable for TCN exposure dating [START_REF] Tylmann | Erratics selection for cosmogenic nuclide exposure dating -an optimization approach[END_REF] and to track the last SIS retreat from its maximum limit.

The ice sheet limit of the LLGM in northern Poland is mostly based on the spatial distribution of end moraines, outlets of tunnel valleys or proximal parts of outwash plains and fans [START_REF] Majdanowski | Rozmieszczenie, gęstość i kierunki rynien jeziornych na Niżu Polskim[END_REF][START_REF] Galon | Extents of the Scandinavian glaciations and of their recession stages on the territory of Poland in the light of analysis of the marginal forms of inland ice[END_REF][START_REF] Marks | Last Glacial Maximum in Poland[END_REF]. This limit defines the "young glacial" landscape located to the north and the "old glacial" landscape located to the south of it (Fig. 1B). According to recent reconstructions, the LLGM was diachronous and it occurred during the Leszno/Brandenburg Phase in the western part and during the Poznań/Frankfurt Phase in the eastern part of the study area [START_REF] Marks | Quaternary glaciations in Poland[END_REF][START_REF] Marks | Timing of the Late Vistulian (Weichselian) glacial phases in Poland[END_REF][START_REF] Marks | Last deglaciation of northern continental Europe[END_REF]; Fig. 1B). It is supported by lithostratigraphic evidences which suggest asynchronous maxima of the last SIS across the study area. In north-central Poland two till layers separated by icewedge casts, are correlated with the glacial phases mentioned above [START_REF] Wysota | Record of the Vistula ice lobe advances in the Late Weichselian glacial sequence in north-central Poland[END_REF][START_REF] Narloch | Sedimentological record of subglacial conditions and ice sheet dynamics of the Vistula Ice Stream (north-central Poland) during the Last Glaciation[END_REF]. However, in the former Vistula Ice Lobe area (Fig. 1B), only one till layer has been identified, and OSL dating of fluvial and glaciofluvial sediments suggests that it is correlated with the Poznań Phase (cf. [START_REF] Wysota | Record of the Vistula ice lobe advances in the Late Weichselian glacial sequence in north-central Poland[END_REF][START_REF] Roman | Ice-flow directions of the last Scandinavian Ice Sheet in central Poland[END_REF]. In contrast, in the western part of the study area, the till layer correlated with the Leszno Phase occurs close to the maximum extent of the last SIS, and no basal till layer correlated with the Poznań Phase has been documented there. Therefore, the ice marginal belt formed during the Poznań Phase in western Poland is commonly interpreted as a result of the ice margin standstill, which occurred during the last SIS recession after the LLGM of the Leszno Phase (e.g. [START_REF] Kasprzak | Differentiation mechanism in the formation of marginal zones of Leszno and Poznań Phases of the Last Glaciation, Great Poland Lowland[END_REF][START_REF] Kozarski | Time and dynamics of the last Scandinavian ice-sheet retreat from northwestern Poland[END_REF][START_REF] Kozarski | Deglaciation of northwestern Poland: environmental conditions and geosystem transformation ~20 ka → 10 ka BP[END_REF]. Limits of the last SIS from [START_REF] Rinterknecht | The last deglaciation of the southeastern sector of the Scandinavian ice sheet[END_REF] and [START_REF] Kozarski | Deglaciation of northwestern Poland: environmental conditions and geosystem transformation ~20 ka → 10 ka BP[END_REF].

Up until recently, the age of the maximum extent of the last SIS in the study area was interpreted mainly based on available radiocarbon ages (Table 1 and references therein).

Published 14 C data significant for the interpretation of the last SIS advance and retreat (locations in Fig. 1A) indicate the possible time window for its maximum extent and the beginning of its retreat. Radiocarbon ages from organic deposits resting below the LLGM till layers have been obtained in the vicinity of the last SIS maximum limit in central Poland [START_REF] Stankowska | Maximum extent of the Vistulian ice sheet in the vicinity of Konin, Poland: a geomorphological, sedimentological and radiometric evidence[END_REF][START_REF] Gogołek | Vistulian of the Turek Upland in the light of new data[END_REF] as well as in the area of the northern Polish coastline [START_REF] Rotnicki | The Last Cold Period in the Gardno-Łeba Coastal Plain[END_REF] and in offshore boreholes drilled into the Baltic Sea bottom [START_REF] Kramarska | Origin and development of the Odra bank in the light of the geologic structure and radiocarbon dating[END_REF][START_REF] Krzyszkowski | Stratigraphy of Weichselian deposits in the cliff sections between Łukęcin and Niechorze, Baltic coast, Northwestern Poland[END_REF]. The youngest calibrated radiocarbon ages of these deposits constrain the maximum reliable age of the LLGM. The most reliable ages for the ice sheet retreat, may be estimated by the oldest calibrated radiocarbon ages of the organic sediments deposited after deglaciation at the bottom of lake basins close to the last SIS maximum extent in central and northern Poland [START_REF] Niewiarowski | Fluctuations of the water level in the Biskupin lake and their reasons[END_REF][START_REF] Stankowski | Upper Quaternary deposits sites Mikorzyn and Sławoszewek in the light of geological investigation, radiocarbon and thermoluminescence dating, and palynological analysis[END_REF][START_REF] Gamrat | Asynchronous development of two Late Glacial lake basins near the Drwęca icemarginal valley (N Poland)[END_REF]. Because of the risk of organic matter redeposition at the bottom of lake basins and the associated overestimated radiocarbon ages, these ages should be interpreted as a maximum reliable time interval for the deglaciation of north-central Poland after the LLGM.

The calibrated radiocarbon ages of organic deposits underlying the LLGM till layers range between 26.5 ± 1.3 and 25.2 ± 3.4 cal ka BP. The oldest minimum age predating the LLGM was found on the northern Polish coastline [START_REF] Rotnicki | The Last Cold Period in the Gardno-Łeba Coastal Plain[END_REF]) and the youngest one in central Poland [START_REF] Gogołek | Vistulian of the Turek Upland in the light of new data[END_REF]. The post-LLGM maximum ages of organic sediments deposited after deglaciation range between 21.4 ± 1.9 cal ka BP and 19.1 ± 0.2 cal ka BP. The oldest age (21.4 ± 1.9 cal ka BP) from bottom part of the lake deposits resting above the LLGM till in central Poland [START_REF] Stankowski | Upper Quaternary deposits sites Mikorzyn and Sławoszewek in the light of geological investigation, radiocarbon and thermoluminescence dating, and palynological analysis[END_REF]. A similar age (21.4 ± 0.6 cal ka BP) was obtained from the bottom of the lacustrine deposits overlying the till of the same stratigraphic position, ~60 km further to the north [START_REF] Niewiarowski | Fluctuations of the water level in the Biskupin lake and their reasons[END_REF].

The youngest maximum post-LLGM age (19.1 ± 0.2 cal ka BP), was found in lacustrine silt filling glacial channel located further to the east [START_REF] Gamrat | Asynchronous development of two Late Glacial lake basins near the Drwęca icemarginal valley (N Poland)[END_REF].

The distribution of pre-LLGM and post-LLGM calibrated 14 C ages in the study area suggests that:

1) the LLGM ice advance occurred after ~25 cal ka BP, 2) the beginning of ice sheet retreat after the LLGM in central-western Poland occurred not earlier than ~21.5 cal ka BP and the deglaciation of the regions located further to the east started not earlier than ~19.0 cal ka BP.

Methods

Sample selection

The largest erratic boulders located in situ on glacial landforms in the vicinity of the maximum extent of the last SIS and to the north of this limit were our targets for sampling for surface exposure dating with TCN (Fig. 1B). The boulders were carefully selected from an extensive dataset of almost 500 erratics, based on their dimensions and geomorphological locations [START_REF] Tylmann | Erratics selection for cosmogenic nuclide exposure dating -an optimization approach[END_REF]. We selected erratics located on depositional glacial landforms, i.e. terminal moraines, dead-ice moraines, moraine plateaux or proximal zones of outwash fans/plains. Samples were collected with a manual jackhammer from the upper surface of 33 stable and massive boulders of perimeter ranging from 6.0 to 14.3 m and of height ranging from 0.8 to 2.9 m (Table 2; see examples in Fig. 2). Almost all boulders (n = 32) are characterized by quartz-rich lithologies -granitoid, granite gneiss and gneiss. Because of the availability of quartz, these erratics were sampled for surface exposure dating with in situ-produced 10 Be. Out of 33 samples for 10 Be analysis (Table 2), one boulder was sampled twice; we took a sample from the top (LGM-02T) and one from the side (LGM-02S). The motivation of this sampling was that from the field observations the boulder looked to have been rotated by 90°, with the original top surface on the side (Fig. 2C). One boulder does not contain quartz (gabbro) and it was sampled for exposure dating with in situ-produced 36 Cl (Fig. 2D, Tables 345). In total, 34 new samples were analyzed: 33 for 10 Be analysis and one for 36 Cl analysis.

10 Be dating methodology

The samples were crushed and sieved. The 0.25-1.0 mm fraction was divided with a Frantz separator into magnetic and non-magnetic subsamples. The non-magnetic fractions were enriched in quartz decontaminated by successive acid leaching (HCl + H 2 SiF 6 then dilute HF). The purified quartz was spiked with ~100-300 mg of a 3025 ppm home-made carrier then dissolved in 48% HF. After evaporation and purification with Dowex 1 X 8 and Dowex 50W * 8 (100-200 mesh) ion exchange resins, BeOH was dried and then oxidized to BeO at 800°C. BeO was mixed with 325 mesh Nb-powder prior to measurement at the accelerator mass spectrometer (AMS). Sample preparation was conducted at the Laboratoire National des Nucléides Cosmogéniques (LN 2 C) at CEREGE, Aix-en-Provence, France. The French AMS facility ASTER was used for the 10 Be measurements [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: Improved performance and developments[END_REF]. The measured 10 Be/ 9 Be ratios were normalized relative to the in-house standard STD-11 using an assigned 10 Be/ 9 Be ratio of (1.191 ± 0.013) × 10 -11 [START_REF] Braucher | Preparation of ASTER in-house 10 Be/ 9 Be standard solutions[END_REF] and a 10 Be half-life of (1.387 ± 0.012) × 10 6 years [START_REF] Chmeleff | Determination of the 10 Be half-life by multicollector ICP-MS and liquid scintillation counting[END_REF][START_REF] Korschinek | A new value for the half-life of 10 Be by Heavy-Ion Elastic Recoil Detection and liquid scintillation counting[END_REF]. 10 Be ages were calculated using the most recent global production rate [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF]) and the 'Lm' time-dependent scaling scheme for spallation according to [START_REF] Lal | Cosmic ray labeling of erosion surfaces: In situ nuclide production rates and erosion models[END_REF] and [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF]. The modern global 10 Be production rate at sea level and high-latitude is derived from a recalibrated global data set and recent calibration studies [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF]. This production rate was used as the most adequate, because the local production rate for Poland was not estimated so far. A recently published Swedish 10 Be production rate [START_REF] Stroeven | A new Scandinavian reference 10 Be production rate[END_REF] would generate slightly younger exposure ages (~5%). The 10 Be production rate was corrected for sample thickness according to an exponential function [START_REF] Lal | Cosmic ray labeling of erosion surfaces: In situ nuclide production rates and erosion models[END_REF] and assuming an average density of 2.7 g/cm 3 . A correction for self-shielding (direction and angle of surface dipping) was applied when the surface of a sampling spot was sloping >10°. However, for most of the sampled boulders no correction for surface inclination was necessary. We apply a correction for surface erosion, assuming a maximum erosion rate of 1.3 mm ka -1 for granitic surfaces [START_REF] Gosse | Precise cosmogenic 10 Be measurements in western North America: support for a global Younger Dryas cooling event[END_REF][START_REF] Rinterknecht | Late Quaternary ice sheet extents in northeastern Germany inferred from surface exposure dating[END_REF][START_REF] Rinterknecht | Cosmogenic 10 Be dating of ice sheet marginal belts in Mecklenburg-Vorpommern, Western Pomerania (northeast Germany)[END_REF]. Corrections for snow/vegetation cover were not included due to their minimal impact on 10 Be ages in the temperate vegetation zone (e.g. [START_REF] Plug | Attenuation of cosmic ray flux in temperate forests[END_REF][START_REF] Rinterknecht | Cosmogenic 10 Be dating of ice sheet marginal belts in Mecklenburg-Vorpommern, Western Pomerania (northeast Germany)[END_REF]. All calculations were performed using the online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator -version 3 (http://hess.ess.washington.edu/math/; accessed: 2019-02-01), which is an updated version of the online calculator described by [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF]. 10 Be exposure ages corrected and uncorrected for erosion are reported with 1σ uncertainties (including analytical uncertainties and the production rate uncertainty) in ka in Table 2. In the following sections we use ages uncorrected for the erosion and we interpret them as minimum ages for the ice sheet retreat.

36 Cl dating methodology

Gabbro sample FRT-11 was crushed and sieved to a grain size fraction of 250-710 µm at CALM (Cosmonucléides Au Laboratoire de Meudon) at the Laboratoire de Géographie Physique (LGP), France, and chemically prepared for in situ 36 Cl dating at the LN 2 C. The procedure is described in [START_REF] Schimmelpfennig | Calibration of cosmogenic 36 Cl production rates from Ca and K spallation in lava flows from Mt. Etna (38°N, Italy) and Payun Matru (36°S, Argentina)[END_REF]. Sample specific data and compositional results are given in Tables 3, 4 and 5. The 36 Cl age of sample FRT-11 was calculated with the Excel® spreadsheet of [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF], using the scaling method of [START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF] and employing the following 36 Cl production rates, referenced to sea level and high latitude (SLHL): 42.2 ± 4.8 atoms 36 Cl (g Ca) -1 yr -1 for spallation of Ca [START_REF] Schimmelpfennig | Calibration of cosmogenic 36 Cl production rates from Ca and K spallation in lava flows from Mt. Etna (38°N, Italy) and Payun Matru (36°S, Argentina)[END_REF], 148.1 ± 7.8 atoms 36 Cl (g K) -1 yr -1 for spallation of K (Schimmelpfennig et al., 2014), 13 ± 3 atoms 36 Cl (g Ti) -1 yr -1 for spallation of Ti [START_REF] Fink | Cross-sections for 36 Cl from Ti at Ep=35-150 MeV: Applications to in-situ exposure dating[END_REF], 1.9 ± 0.2 atoms 36 Cl (g Fe) -1 yr -1 for spallation of Fe [START_REF] Stone | Terrestrial chlorine-36 production from spallation of iron[END_REF], and 696 ± 185 neutrons (g air) -1 yr -1 for the production rate of epithermal neutrons from fast neutrons in the atmosphere at the land/atmosphere interface (Marrero et al., 2016b). We used a high-energy neutron attenuation length of 160 g cm -2 . The resulting 36 Cl age in Table 5 with its 1 σ uncertainty including all analytical and production rate errors is shown with and without correction for an erosion rate of 1.3 mm ka -1 . The 36 Cl age was also calculated using the CRONUScalc calculator (Marrero et al., 2016a) for comparison with the results calculated with the Excel® spreadsheet [START_REF] Schimmelpfennig | Sources of in-situ 36 Cl in basaltic rocks. Implications for calibration of production rates[END_REF].

Results

Erratics located beyond the ice marginal belt of the LLGM

Exposure ages for SAAL samples show significant scatter (16.2 ± 1.4 to 115.2 ± 11.0 ka; Table 2). The observed variability (53.2%) is much larger than the average external uncertainty (8.4%). This shows that random uncertainties are dominated by geological uncertainties rather than by analytical ones. Most of the pre-LLGM exposure ages are younger than 100 ka (six samples) with five ages younger than 60 ka (Table 2).

Erratics of the LLGM ice marginal belt

We sampled 13 erratics located on the ice marginal belt identified geomorphologically as the Leszno (Brandenburg) Phase (Fig. 1B). The exposure ages for samples LGM (01 to 08) and FRT (01 to 04 and 06) range between 19.3 ± 2.1 and 75.3 ± 6.3 ka (Table 2). The distribution of the ages is bimodal with the narrow and high mode containing ages ranging between 19.3 ± 2.1 and 25.1 ± 2.1 ka, and wider and lower mode containing ages between 39.5 ± 3.3 and 75.3 ± 6.3 ka (Fig. 3A). The results for the boulder that has been sampled twice (LGM-02) belong to the older distribution, and are: 39.5 ± 3.3 ka for the top sample (LGM-02T) and 58.6 ± 3.9 ka for the side sample (LGM-02S). The older age of the weathered side of the boulder support our hypothesis that this erratic was rotated and the original top surface is on the side (Fig. 2C). Ages belonging to the younger distribution are much more clustered around the mode than the ages of the older distribution (Fig. 3A). The reduced chi-squared test for ages of the younger distribution indicates that they represent a single population (χ 2 r = 0.8 and K = 2.2), while the ages of the older distribution do not represent a single population (χ 2 r = 7.2 and K = 2.2). We used the criterion 'K' to test whether the reduced chi-squared value fall within a 2σ envelope (95% confidence). This criterion depends on the degree of freedom, and thus the number of samples (n): K = 1 + 2√2/(n -1). If χ 2 r < K then there is a >95% probability that the analyzed measurements represent a single population. We therefore used the younger distribution to interpret the age of the Leszno Phase ice marginal belt. The oldest age of this dataset (25.1 ± 2.1 ka) was identified as an outlier according to the Chauvenet`s criterion. A Shapiro-Wilk test, applied after rejecting the outlier, showed that we cannot reject the normality assumption for p = 0.05 (W = 0.97) for distribution of the remaining six ages. The error-weighted mean and the standard deviation of the mean calculated for these ages is 20.7 ± 0.8 ka (Table 2). The spatial distribution of the analyzed erratics indicates that statistically reliable ages (included in the mean age calculation) are located mainly in the western part of the study area (Fig. 4).

Twelve boulders located on the Poznań (Frankfurt) ice marginal belt were sampled for 10 Be (samples LGM-09 to12, FRT-05, FRT-07 to 10 and FRT-12 to 14) and one boulder for 36 Cl exposure dating (sample FRT-11). We also recalculated and used 10 Be ages of ten boulders from published data [START_REF] Rinterknecht | Cosmogenic 10 Be ages on the Pomeranian moraine[END_REF][START_REF] Rinterknecht | The last deglaciation of the southeastern sector of the Scandinavian ice sheet[END_REF], located close to the ice sheet limit during the LLGM in the eastern edge of the study area (Table 2; Fig. 1B) using the same procedure as for our samples. 10 Be ages range between 5.8 ± 0.8 and 40.6 ± 3.9 ka (Table 2) with narrow modes occurring at ~5-7 ka and ~15-19 ka, and wide, low mode at ~30-40 ka (Fig. 3B). The youngest mode includes two ages 5.8 ± 0.8 and 7.0 ± 1.9 ka, the intermediate mode consists of 14 ages ranging between 15.5 ± 1.6 and 20.2 ± 2.2 ka, and the oldest mode includes four ages between 30.2 ± 3.2 and 40.6 ± 3.9 ka. The reduced chi-squared test indicates that the 22 10 Be ages represent multiple populations (χ 2 r = 18.0 and K = 1.6). The application of Chauvenet`s criterion identifies that the three youngest ages (<13.0 ka) and the five oldest ages (>24.1 ka) are outliers (Fig. 3B). After rejecting the outliers, the reduced chi-squared test applied to the remaining 14 ages indicates that they belong to a single population (χ 2 r = 0.5 and K = 1.8). Moreover, a Shapiro-Wilk test showed that we cannot reject the normality assumption for p = 0.05 (W = 0.90) for the distribution of the 14 ages. The error-weighted mean and the standard deviation of the error-weighted mean calculated for the 10 Be ages is 17.3 ± 0.5 ka (Table 2). The 36 belong to the same frequency distribution as the 14 10 Be ages (Fig. 3B).

Fig. 4. Spatial distribution of 10 Be ages and single 36 Cl age in the study area. Ages used in the calculation of the error-weighted mean ages for specific ice marginal belts are framed in black frame.

Discussion

Pre-LLGM erratics

The distribution of the exposure ages of erratics located beyond of the ice sheet limit during the LLGM is very large (from 15.9 ± 1.4 to 101.7 ± 8.5 ka; Table 2) and these ages are not clustered around any particular time interval. We have only one exposure age (sample SAAL-05) that may be roughly interpreted as representing the timing of retreat of the MIS 6 ice sheet: 101.7 ± 8.5 ka (Table 2). Six out of seven boulders display much too young ages (from 15.9 ± 1.4 to 66. Two erratics (samples SAAL-04: 22.3 ± 1.9 ka and SAAL-06: 15.9 ± 1.4 ka) have exposure ages close to the expected LLGM age, and they are also located close to the ice sheet limit during the LLGM (Figs. 1 and4). This may suggest a revision of the last SIS maximum extent in central Poland. However, these erratics are located in the vicinity of the sites of Eemian organic deposits that are not covered by the LLGM till [START_REF] Tobolski | Paleobotanical studies of the Eemian interglacial and early Vistulian, Władysławów in the vicinity of Turek (preliminary report)[END_REF][START_REF] Kuszell | New sites of the interglacial deposits in south-western Poland[END_REF]; Szałamacha and Skompski, 1999; Balwierz and Roman, 2002; Fig. 1). In addition to geomorphology, this is a strong argument supporting that the ice sheet limit during the LLGM is correct and that these two boulders display underestimated exposure ages resulting from geomorphic processes.

Age of the LLGM and the first stages of the last SIS retreat

In central-western and northern Poland, deposition of organic sediments resting below the LLGM till is dated to ~26.5-25.2 cal ka BP (Table 1) indicating that the maximum extent of the last SIS in this region occurred later than ~25 cal ka BP. Moreover, the radiocarbon age of the basal part of gyttja (21.4 ± 1.9 cal ka BP) resting on the same till [START_REF] Stankowski | Upper Quaternary deposits sites Mikorzyn and Sławoszewek in the light of geological investigation, radiocarbon and thermoluminescence dating, and palynological analysis[END_REF] indicates that the ice sheet retreat from its maximum limit started not earlier than ~21.4 cal ka BP. Recently, Marks (2012) determined, based on his interpretation of the available ages, the age of the maximum extent of the last SIS in central and western Poland during the Leszno Phase at 24 cal ka BP.

Our results for the Leszno (Brandenburg) Phase in western Poland are clustered around two modes at ~19-22 ka and ~40-70 ka (Fig. 3A). We used this younger mode to calculate an error-weighted mean surface exposure age of 20.7 ± 0.8 ka (Fig. 5). We interpret this age as a minimum age of the commencement of the ice sheet retreat from its maximum limit at the end of the LLGM in central and western Poland. Ages used in the error-weighted mean age calculation agree, within uncertainty, with the radiocarbon ages coming from deposits resting below and above the till of the ice advance during the Leszno Phase in central Poland [START_REF] Stankowski | Upper Quaternary deposits sites Mikorzyn and Sławoszewek in the light of geological investigation, radiocarbon and thermoluminescence dating, and palynological analysis[END_REF]. The ages of the older distribution (between 39.5 ± 3.3 and 75.3 ± 6.3 ka) probably result from inherited 10 Be. This could have happened either due to redeposition of previously exposed boulders by the last SIS with no or limited erosion (cf. [START_REF] Håkansson | Late Pleistocene glacial history of Jameson land, central east Greenland, derived from cosmogenic 10 Be and 26 Al exposure dating[END_REF][START_REF] Rinterknecht | Late Quaternary ice sheet extents in northeastern Germany inferred from surface exposure dating[END_REF].

The new exposure age of the last ice sheet retreat in western Poland is also consistent with most luminescence ages constraining the timing of the Leszno Phase (Fig. 5). OSL Poland, constrained with the available radiocarbon ages and our new 10 Be exposure ages, is between ~25 ka and ~21 ka (Fig. 5).

Radiocarbon dating from western Belarus [START_REF] Zimenkov | Age of the maximum extent of the Poozerian Glaciation in Belarus. New data about Cenozoic geology of Belarus and adjacent areas[END_REF] indicates an age of the last SIS maximum extent of ~22.5 cal ka BP in the eastern part of its southern sector. The deglaciation of this region, in the light of radiocarbon dating in north-eastern Poland [START_REF] Gamrat | Asynchronous development of two Late Glacial lake basins near the Drwęca icemarginal valley (N Poland)[END_REF], started not earlier than 19.1 ± 0.2 cal ka BP. Our results for erratics located close to the ice sheet limit during the LLGM in the eastern part of the study area also support this view (Fig. 5). TCN exposure ages of boulders located on the ice marginal belt of the Poznań Phase, after rejecting outliers with Chauvenet`s criterion, range between 15.5 ± 1.6 and 20.2 ± 2.2 ka, with an average age of 17.3 ± 0.5 ka (Table 2). We interpret these ages as a minimum age of: ± 0.8 and 20.1 ± 1.2 ka respectively. Furthermore, the recent results from Valday Heights in NW Russia show a 10 Be age for the LLGM at 18.9 ± 0.5 ka [START_REF] Rinterknecht | The Last Glacial Maximum extent of the Scandinavian Ice Sheet in the Valday Heights, western Russia: Evidence from cosmogenic surface exposure dating using 10 Be[END_REF]. These data suggest an older age for the LLGM further east of eastern Poland. We propose a possible duration for the LLGM in eastern Poland, constrained by radiocarbon ages and our new 10 Be exposure ages as well as recalculated data of [START_REF] Rinterknecht | Cosmogenic 10 Be ages on the Pomeranian moraine[END_REF][START_REF] Rinterknecht | The last deglaciation of the southeastern sector of the Scandinavian ice sheet[END_REF], between ~22 ka and ~18 ka (Fig. 5). Table 1. Published radiocarbon ages significant for the timing of the last SIS advance and retreat along the LLGM limit in Poland.

Fig. 1 .

 1 Fig. 1. The main ice margin limits of the last Scandinavian Ice Sheet (SIS) in the area south of

Fig. 2 .

 2 Fig. 2. Examples of erratic boulders sampled for TCN exposure dating. (A) Sampling of FRT-

Fig. 3 .

 3 Fig. 3. Distribution of surface exposure ages for erratic boulders located on the LLGM ice

  Cl age is 20.1 ± 3.1 ka when calculated with the Excel® spreadsheet (Schimmelpfennig et al., 2009), and 20.0 ± 3.0 ka when calculated with the CRONUScalc calculator (Marrero et al., 2016b). These two results are very similar and

  8 ± 5.5 ka) for the retreat of the MIS 6 ice sheet, and probably result from degradation of moraines and erosion of erratics (e.g.[START_REF] Hallet | Surface exposure dating of dynamic landforms: young boulders on aging moraines[END_REF] during the long time period after deglaciation. Postglacial weathering and erosion of an erratic surface, removes cosmogenic nuclide accumulated in the surface layer of the sample and thus produces too young ages, but we consider the obtained exposure ages most probably represent long-term geomorphic processes acting after deglaciation, and therefore largely underestimate the timing of ice sheet retreat. Similar inconsistencies of10 Be ages for erratics located within Saalian (MIS 6) landscapes were reported in Germany, in the Hocher Fläming area[START_REF] Rinterknecht | Late Quaternary ice sheet extents in northeastern Germany inferred from surface exposure dating[END_REF], and in Mecklenburg-Vorpommern[START_REF] Rinterknecht | Cosmogenic 10 Be dating of ice sheet marginal belts in Mecklenburg-Vorpommern, Western Pomerania (northeast Germany)[END_REF] in front of the ice sheet limit during the Brandenburg Phase. For 14 erratics sampled for10 Be dating in Germany, ten boulders have exposure ages too young (from 27.7 ± 1.3 to 92.8 ± 3.0 ka) for the deglaciation at the end of MIS 6 and four erratics have exposure ages ranging from 101.7 ± 10.4 to 143 ± 5.0 ka.

  dating of fluvial, fluvioglacial and periglacial sediments underlying the till correlated with the Leszno Phase (Stankowski et al., 1999; Wysota et al., 2009; Molewski and Weckwerth, 2009; Weckwerth et al., 2011) yields ages that mostly range between ~29 and ~24 ka. Additionally, OSL ages of sediments underlying the LLGM till in the northernmost region of Germany suggest that the last ice sheet reached the southern Baltic coast at Rügen Island after 22.7 ± 1.9 ka (Pisarska-Jamroży et al., 2018) or 22.0 ± 2.0 ka and 23.0 ± 2.0 ka (Kenzler et al., 2015, 2017). OSL ages of fluvial and periglacial deposits overlying the till correlated with the Leszno Phase in Poland range between ~21 and 19 ka (Molewski 2007; Wysota et al., 2009; see Fig. 5). We thus suggest that the error-weighted mean exposure age of erratics in western Poland (20.7 ± 0.8 ka) represents the most reliable age for the initiation of the ice sheet retreat after the ice advance during the Leszno Phase. Our results also corroborate TCN data published by Heine et al. (2009), who determined the exposure 10 Be age of three boulders located on the Brandenburg (Leszno) and Frankfurt (Poznań) Phase ice marginal belts in north-eastern Germany. Our error-weighted mean exposure age (20.7 ± 0.8 ka) and the results of Heine et al. (2009) (22.0 ± 1.1 ka), are coeval within uncertainties. We suggest that the duration of the LLGM in central and western

1 )

 1 the beginning of the ice sheet retreat during the LLGM in eastern Poland, after the ice advance during the Poznań Phase, 2) the ice sheet retreat from the Poznań Phase limit in north-central Poland, postdating the LLGM represented by the Leszno Phase. These results are in general agreement with OSL ages (~18.7-16.7 ka; Fig. 5) of fluvioglacial sand overlying the till correlated with the ice advance during the Poznań Phase in north-central Poland[START_REF] Stankowski | Upper Quaternary deposits sites Mikorzyn and Sławoszewek in the light of geological investigation, radiocarbon and thermoluminescence dating, and palynological analysis[END_REF][START_REF] Wysota | Record of the Vistula ice lobe advances in the Late Weichselian glacial sequence in north-central Poland[END_REF][START_REF] Roman | Ice-flow directions of the last Scandinavian Ice Sheet in central Poland[END_REF].10 Be exposure dating of erratics located on moraines correlated with the LLGM in Belarus and Lithuania[START_REF] Rinterknecht | Timing of the Last Deglaciation in Belarus[END_REF][START_REF] Rinterknecht | Timing of the last deglaciation in Lithuania[END_REF] yields mean exposure ages of deglaciation at 17.7

Fig. 5 . 5 . 3 . 7 ± 0.8 ka. 3 . 6 .

 553736 Fig. 5. Distribution of TCN, 14 C, and OSL ages relevant for the reconstruction of the timing of the LLGM and the initiation of the last ice sheet

Table 2 .

 2 Surface exposure10 Be ages of pre-LLGM and LLGM erratic boulders in northern Poland. The list consists of 33 new exposure ages and 10 10 Be ages recalculated from the original data of[START_REF] Rinterknecht | Cosmogenic 10 Be ages on the Pomeranian moraine[END_REF][START_REF] Rinterknecht | The last deglaciation of the southeastern sector of the Scandinavian ice sheet[END_REF]. Densities of all samples (granitoids, granite gneisses and gneisses) were assumed at 2.7 g/cm 3 . All 10 Be exposure ages are calculated with 'Lm' time-dependent scaling scheme for spallation according to[START_REF] Lal | Cosmic ray labeling of erosion surfaces: In situ nuclide production rates and erosion models[END_REF] and[START_REF] Stone | Air pressure and cosmogenic isotope production[END_REF] and global production rate according to[START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF].

	Laboratory symbol 14 C age (ka BP)	Calibrated	Dated material	Locality	Reference
			age (cal ka
			BP)*
	Post-LLGM ages
	Poz-46881	15.9 ± 0.1	19.1 ± 0.2	dispersed organic matter in laminated lacustrine silt	Piotrowo (Poland)	Gamrat et al. (2017)
	Gd-5777	17.7 ± 0.2	21.4 ± 0.6	dispersed organic matter in lacustrine silt	Biskupińskie Lake (Poland)	Niewiarowski (1995)
	Gd-9604	17.7 ± 0.8	21.4 ± 1.9	peat at the bottom of gyttja	Sławoszewek (Poland)	Stankowski et al. (1999)
	Pre-LLGM ages
	Gd-786	21.2 ± 1.6	25.2 ± 3.4	dispersed organic matter in silt below till	Żdzary (Poland)	Gogołek and Mańkowska (1989)
	Gd-1038	21.5 ± 0.4	25.6 ± 1.1	dispersed organic matter in silt below till	Odra Bank (Poland)	Kramarska (1998)
	Gd-9686	21.6 ± 1.1	25.6 ± 2.2	organic matter below till	Niechorze (Poland)	Krzyszkowski et al. (1999)
	Gd-645	22.1 ± 0.5	26.4 ± 0.9	peat below till	Maliniec (Poland)	Stankowska and Stankowski (1988)
	Gd-4581	22.3 ± 0.7	26.5 ± 1.3	dispersed organic matter in silt	Łeba Barrier (Poland)	Rotnicki and Borówka (1995)

* Calibrated using IntCal13 calibration curve (Reimer et al., 2013) with the Calib 7.0.4 software (Stuiver et al., 2017).

.9^ Poznań (Frankfurt) Phase Mean 10 Be age (n = 14): 17.3 ± 0.5

	LGM-09	52.7390	19.6783	143	14.3/1.4	granitoid	1.0	1.0000	15.21	0.1036	15.93 ± 0.77	34.3 ± 3.2	33.0 ± 3.0
	LGM-10	52.8749	19.6290	114	9.0/0.8	granitoid	1.0	1.0000	17.31	0.1021	7.95 ± 0.38	17.3 ± 1.6	17.0 ± 1.5L
	GM-11	53.4187	19.9013	166	9.5/1.2	granitoid	1.1	1.0000	14.68	0.1032	8.67 ± 0.47	17.8 ± 1.7	17.5 ± 1.6L
	GM-12	53.3874	19.5752	130	8.0/1.0	granitoid	1.2	1.0000	15.07	0.1026	11.50 ± 0.53	24.8 ± 2.3	24.1 ± 2.1
	FRT-05	52.4897	15.1843	113	8.7/1.6	granitoid	1.1	1.0000	15.60	0.1023	14.13 ± 1.06	31.3 ± 3.5	30.2 ± 3.2
	FRT-07	53.2823	17.5074	116	8.9/1.2	granitoid	1.0	1.0000	9.71	0.3080	8.91 ± 1.07	19.3 ± 2.8	18.9 ± 2.7F
	RT-08	52.9537	17.1255	89	17.5/2.2	gneiss	2.5	1.0000	15.03	0.1023	7.74 ± 0.43	17.5 ± 1.7	17.1 ± 1.6F
	RT-09	52.9525	17.3370	94	13.5/1.3	granitoid	1.2	1.0000	15.65	0.1025	7.66 ± 0.38	17.0 ± 1.6	16.7 ± 1.5F
	RT-10	52.7259	17.3231	105	20.5/1.1	granitoid	0.7	1.0000	16.46	0.1025	7.89 ± 0.65	17.3 ± 2.0	16.9 ± 1.9F
	RT-12	53.5238	16.2029	150	19.0/3.4	gneiss	1.6	0.9996	22.03	0.3049	9.81 ± 0.76	20.6 ± 2.3	20.2 ± 2.2F
	RT-13	53.6903	17.1372	171	20.5/1.9	gneiss	2.0	1.0000	9.29	0.3055	9.69 ± 0.57	20.0 ± 2.0	19.5 ± 1.9F
	RT-14	53.7271	17.1342	187	14.6/2.1	granitoid	2.5	0.9994	20.20	0.3052	8.55 ± 0.54	17.4 ± 1.7	17.0 ± 1.7L
	ES-3	53.4806	20.2222	212	n.a./1.6	granitoid	2.0	1.0000	39.99	0.2500	3.45 ± 0.35	7.1 ± 0.9	7.0 ± 0.9
	LES-4	53.4944	20.2389	172	n.a./1.8	granitoid	2.0	1.0000	40.00	0.2500	7.40 ± 0.54	15.8 ± 1.7	15.5 ± 1.6*L
	ES-5	53.5792	20.0944	180	11.7/2.5	granitoid	2.0	1.0000	40.00	0.2500	19.24 ± 1.16	42.5 ± 4.3	40.6 ± 3.9*
	LES-6	53.6006	20.0611	151	17.1/2.2	gneiss	2.0	1.0000	40.00	0.2500	8.08 ± 0.58	17.7 ± 1.9	17.4 ± 1.8*L
	ES-7	53.6250	20.0417	132	13.0/2.8	granitoid	2.0	1.0000	60.51	0.2500	2.64 ± 0.33	5.8 ± 0.9	5.8 ± 0.8*
	LES-8	53.5111	19.9000	255	4.5/1.1	granitoid	2.0	1.0000	40.00	0.2500	10.14 ± 1.10	20.1 ± 2.7	19.7 ± 2.6*L
	ES-10	53.5764	19.9417	270	8.7/1.7	granitoide	2.0	1.0000	40.01	0.2500	6.78 ± 0.57	13.2 ± 1.5	13.0 ± 1.5*
	LES-11	53.5222	19.8375	218	10.5/1.5	gneiss	2.0	1.0000	39.99	0.2500	7.94 ± 0.77	16.3 ± 2.0	16.0 ± 2.0*L
	ES-12	53.5625	19.9528	275	10.1/1.6	granitoid	2.0	1.0000	40.01	0.2500	8.46 ± 0.70	16.4 ± 1.9	16.1 ± 1.8*L
	ES-13	53.5528	19.9250	302	11.0/1.2	granitoid	2.0	1.0000	40.01	0.2500	19.15 ± 1.33	37.0 ± 4.0	35.5 ± 3.7*

Acknowledgements

We are very grateful to the Regional Directorates of Environmental Protection local communes offices for permissions for sampling large erratics protected by law. We appreciate suggestions of two anonymous reviews and of the editor which help to improve the manuscript.

Research is funded by:

National Science Centre grant no. 2014/15/D/ST10/04113 to Karol Tylmann. The ASTER AMS national facility (CEREGE, Aix en Provence) is supported by the INSU/CNRS, the ANR through the "Projets thématiques d'excellence" program for the "Equipements d'excellence" ASTER-CEREGE action and IRD.

All samples measured at the ASTER facility. AMS results are standardized to STD-11 (samples LGM and FRT) and SRM 4325 (samples LES). 10 Be/ 9 Be ratios were corrected for a process blank value of 6.08 x 10 -12 (samples LGM and FRT-01 to FRT-10), 5.97 x 10 -12 (samples FRT-11 to FRT-14), and 26.8 ± 1.4 x 10 -12 (samples LES). 1 Corresponding to self-shielding (direction and angle of surface dipping). * Surface exposure ages recalculated from the original data of [START_REF] Rinterknecht | Cosmogenic 10 Be ages on the Pomeranian moraine[END_REF][START_REF] Rinterknecht | The last deglaciation of the southeastern sector of the Scandinavian ice sheet[END_REF]. ^ Ages used in the error-weighted mean age calculation for the LLGM moraines. n.a. -not available Table 3. Bulk composition of sample FRT-11 before chemical treatment, analysed at the SARM-CRPG (Nancy, France) by ICP-OES (major elements), ICP-MS (trace element), atomic absorption (Li), colorimetry (B) and spectrophotometry (Cl). FRT-11 50. [START_REF] Mix | Environmental processes of the ice age: land, oceans, glaciers (EPILOG)[END_REF]

Sample Name