{'fraction_non_alphanumeric': 0.03875968992248062, 'fraction_numerical': 0.019564414913252122, 'mean_word_length': 4.454697986577181, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalable Parallel Numerical CSP Solver

We present a parallel solver for numerical constraint satisfaction problems (NCSPs) that can scale on a number of cores. Our proposed method runs worker solvers on the available cores and simultaneously the workers cooperate for the search space distribution and balancing. In the experiments, we attained up to 119-fold speedup using 256 cores of a parallel computer.


Introduction
Numerical constraint satisfaction problems (NCSPs, Section 2) have been successfully applied to problems described in the domain of reals [3,6]. Given a NCSP with search space represented as a box (i.e., interval vector), the branch and prune algorithm efficiently computes a paving, a set of boxes that encloses the solution set, yet its exponential computational complexity limits the tractable instances. Although the solving process exhibits a parallelism, no parallel NCSP solver has been made available to date because of the difficulty in partitioning the search space equally.
In this research, we parallelize a NCSP solver to scale its solving process on both shared memory and distributed memory parallel computers (Section 3). Our parallel method consists of parallel worker solvers that solve a portion of search space on CPU cores and interact with neighbor workers via message passing for dynamic load balancing. We also propose a preprocess that accelerates the initial search space distribution by sending sets of boxes via static routing between the workers. We have implemented the method by extending the Realpaver solver using the X10 language to realize a process-level parallelization over a number of cores. Section 4 reports experimental results when our method was deployed on two hardware environments.
Related work. There have been several work regarding parallel solving of CSPs with either discrete or continuous domains. Parallel solving of generic CSPs on massive computer clusters and supercomputers has been explored in [12,2,18]. This work focuses on a massive parallel solver for NCSPs that has a different characteristics compared to generic CSPs. In the survey [7], existing work is classified into (i) search-space splitting methods [16,12,2,5,15,18], (ii) cooperative methods for heterogeneous solver agents (cf. portfolios) [2], and (iii) parallelization of constraint propagation. Our method belongs to the first category. A few works have used approaches (ii) [11] and (iii) [8] for parallelization of NCSP solving. However, to the best of our knowledge, a massive parallelization method that uses the typical approach (i) has not yet been proposed.
Substantial work regarding the parallelization of the branch and bound algorithm with search-space splitting exists [9,13,14]. This approach has also been applied to CSP solvers [16,5] and SAT solvers [15]. This work explores an efficient parallel method for solving NCSPs with similar approach to [9,13].

Numerical Constraint Satisfaction Problems
A numerical constraint satisfaction problem (NCSP) is defined as a triple (v, v 0 , c) that consists of a vector of variables v = (v 1 , . . . , v n ), an initial domain in the form of a box v 0 ∈ IR n (IR denotes the set of closed real intervals), and a con- ≥ 0, where f : R n → R e and g : R n → R i , i.e., a conjunction of e equations and i inequalities. A solution of a NCSP is an assignment of its variables v ∈ v 0 that satisfies its constraints. The solution set Σ of a NCSP is the region within its initial domain that satisfies its constraints, i.e., Σ(v 0 ) := {ṽ ∈ v 0 | c(ṽ)}. The target of this paper is under-constrained NCSPs such that n > e. In general, an under-constrained problem may have a continuous set of infinitely many solutions.
Branch and Prune Algorithm. The branch and prune algorithm [17] is the standard solving method for NCSPs. It takes a NCSP and a precision ǫ as an input and outputs a set of boxes (or paving) S that approximates the solution set with precision ǫ. Examples of S are illustrated in Figure 1.
An intermediate state of the algorithm is represented as a pair of sets of boxes (L, S). The solver receives an initial state ({v 0 }, ∅) and iteratively applies the step computation (illustrated in Figure 2) until it reaches a final state (∅, S ). In the step computation, first, it takes the first element of the queue L of boxes and applies the Prune procedure, which is a filtering procedure that shaves boundary portions of the considered box. In this work, we use an implementation proposed in [6] which provides a verification process based on an interval Newton method combined with a relatively simple filtering process based on the Hull consistency [1]. As a result, a box becomes either empty, precise enough (its width is smaller than ǫ), verified as an inner box of the solution set Σ, or undecided. Precise and inner boxes are appended to S and undecided boxes are passed to Branch. Second, the Branch procedure bisects the box at the midpoint along a component corresponding to one of the variables and the sub-boxes will be put back in the queue. In this work, we assume Branch selects variables in an order which makes the search to behave in a breadth-first manner and thus the solving process gradually improves the precision of the overall pavings ( Figure 1).
The computation of Prune is expensive and is the bottleneck of the solving process. Under certain conditions, application of Prune contracts a large portion of the search space into a tight box (cf. quadratic convergence of the interval Newton methods). Prune can also filter out the whole box if the considered box is verified as an inner or totally inconsistent region. These characteristics of Prune result in the unbalanced nature of search trees. Therefore, a straightforward parallel method does not work efficiently. It is crucial for efficient NCSP solving to execute Prune on each step of traversing the search tree which makes it more difficult to distribute a search path among processors. These properties will be discussed in Section 4.1.

Parallel Branch and Prune
We propose a parallel method that runs several workers on the available CPU cores p 0 , . . . , p #p−1 (we assume a single worker runs on a core and identify both a worker and a core with p i ). Workers homogeneously interleave the following three procedures and cooperate in a decentralized fashion: (i) breadth-first branch and prune search, (ii) distribution and workload balancing of search space in a sender-initiated manner, and (iii) termination detection. Sub-trees in the search space of the branch and prune becomes unbalanced but can be searched independently: there are no confluence of multiple branches, and we have no shared information between branches. Distribution of the search space among workers is done in preprocess and postprocess as described in the following subsections. The preprocess distributes the portions of the search space to the other workers, and the postprocess balances the load of each worker during search. Termination is detected by circulating a termination token via idling workers based on Dijkstra's method (see Section 11.4.4 of [9]).
Preprocess: Search-Space Splitting and Distribution. A solving process is started by a worker p 0 that possesses the initial domain (i.e., a box which contains the whole search space) in the queue L. To distribute the subsequent search space (i.e., a queue of boxes) equally to each core, the preprocess invokes a partition of the queue in two (or more) and then sends a portion to another worker. Figure 3 illustrates in a downward direction some initial transitions of the box queue L distributed among three workers. In our implementation, the distribution routing is formed as a binary tree whose height is ⌈log 2 #p⌉. In each node of the tree, the branch and prune process runs until the number of boxes in the queue reaches nb. Next, the queue is sorted by the volume of the boxes and the half of the content (i.e., nb/2 boxes) is sent to the other core via the right branch.
Postprocess: Dynamic Load Balancing. During search, each worker normalizes the loads within a predefined neighborhood which consists of a small number of neighbor workers. Because there are sufficiently large number of boxes, we simply regard the number of boxes in the queue L as the amount of load. Assume #p workers are running and each worker p i possesses l i boxes in its queue. We also assume for each worker p i that N i is a set of |N | neighbor workers, N −1 i is a set of workers where p j ∈ N −1 i ⇔ p i ∈ N j , L i is a set of loads of the neighbor workers, and ∆ is a predefined load margin. The load balancing procedure of a worker p i performs the following steps once every ns branch and prune steps.
1. For each worker p j in N −1 i , inform the load l i and put in the list L j . 2. Calculate the mean µ of the loads in L i . 3. If µ < ∆, for each worker p j in N i , send at most µ − l j boxes to p j (to be efficient, a certain number of boxes should be kept locally).
Neighbor workers can be identified e.g. as adjacent nodes in the |N |-dimensional mesh of workers. The routing between neighbor workers is fixed during a solving process and thus it may happen that a worker possesses an excessive load than others. However, this load imbalance will be resolved by the subsequent load balancing processes.

Experimental Results
We have implemented the proposed method and measured the speedup of the solving process of under-constrained NCSPs. The experiments were performed with an exhaustive set of parameter combinations to explore the optimal settings. Implementation. We have implemented the proposed method with C++ (gcc ver. 4.4.7 and 4.3.4) and X10 (ver. 2.3.1) [4], a high productivity language for parallel computing. In the following, we use the term place, which is a notion of X10 that in our setting represents a CPU core. Libraries Realpaver (ver. 1.1) [10] for sequential NCSP solving and Gaol (ver. 4.0.1) 4 for basic interval computation were used to facilitate the implementation. The Prune procedure was realized by calling the sequential implementation in Realpaver. Each run of Prune took around 0.2-1ms and the overall execution became the bottleneck of the branch and prune algorithm (occupied greater than 95% of running time in sequential solving). The procedures for search space distribution and load balancing were implemented with X10. Communication of boxes and loads between places were implemented as async tasks and performed in parallel to the search process so that the overhead will be hidden. In the experiments, timings t 1 for sequential runs on single core were measured using the C++ implementation described in [6], which worked identically and faster than our X10 version.
Experiment Environments. Two sets of experiments were operated using (1) a shared-memory machine equipped with 40 cores (four of 10-core Intel Xeon E7-4860 2.26GHz) and 256GB of local memory and (2) up to 256 cores of SGI UV1000, a pseudo-shared memory machine equipped with 2,048 cores (8-core Xeon E7-8837 2.67GHz) and 16TB of memory. UV1000 works as a single shared memory machine by emulating memory accesses using communication based on a high speed NUMAlink5 network which has a bandwidth of 120Gbps. We used the MPI backend of X10 with options X10 NTHREADS = 6 and GC NPROCS = 2.
Experiments on a Shared Memory Machine. We solved the problems shown in [6,3] using 40 cores of the machine (1). We report the results for two representative problems. Parameters in the load balancing method were set as either combination of the following values: nb = 32, |N | ∈ {2, 4}, ns ∈ {10, 100, 1000}, and ∆ = 10. We also computed with and without the preprocess (when the preprocess is not used, the postprocess is executed from the beginning). For each problem, we solved two instances with two multiplicative precisions. The specification of each instance and the computational results are presented in Table 1.
In the table, the columns "problem", "size", "ǫ", "pp", and "|N |" represent the name of the problem, size (i.e., the number of projection/parameter variables), the precision, usage of the preprocess, and the number of neighbor workers, respectively. The rest of the columns represents the results. t 1 and #br 1 represent the running time and the number of branches on single core. t ns,i and #br ns ,i represent the running time and the largest number of branches performed by a worker when computed with the interval ns and i X10 places (best timings are underlined). Figure 4 illustrates the speedup of the solving process.
Experiments on a Cluster with High-speed Interconnection. We solved the problem "3rpr " using up to 256 cores of the machine (2), UV1000. Parameters in the load balancing method were set as either combination of the following values: nb = 8, |N | ∈ {2, 4}, ns = 1000, and ∆ = 10. The results are presented in Table 2. Each column of the table represents the same information as presented in Table 1 except that the column "#sends 1000,256 " represents the number of loads sent by the load balancer in the solving process with ns = 1000 and 256 X10 places. Figure 5 illustrates the speedup of the solving process.

Discussions
In the experiments, our method scaled up to 256 cores with the optimal configurations. We achieved speedups up to 32.3 fold using 40 cores of the shared memory machine and up to 119 fold using 256 cores of the cluster machine. The best speedup of 119 fold was obtained with the preprocess. The preprocess facilitates and accelerates the workload distribution in the early stage of the search process. In some of the experiments without using the preprocess, the speedup ratio became saturated when using many cores (e.g., sp2-4 with ns = 1000 and the experiments on the cluster). This was because the load balancing process was too infrequent for the given number of workers and the work load diffusion became too slow. When comparing the right-hand graphs for the instance sp2-4, ns = 1000, in Figure 4, we can notice that the point of saturation shifts according to the search space size. On the other hand, in some other experiments, the results got worse with the preprocess (e.g., results with ns = 10 on the machine (1)). It occasionally happens that the preprocess mostly solves the problem. However, the preprocess can result in highly unbalanced search trees because of the Prune process, and in such cases the postprocess will not have enough time for load balancing. Regarding the neighborhood sizes |N | = 2, 4, there was a trade off between the workloads balance and the amount of communications required. For the   shared memory machine, it was unclear which size had the advantage. However, for UV1000, the solver was notably slower for |N | = 4 than |N | = 2. It is understandable because larger number of neighbors significantly increased the number of communications (see "#sends 1000,256 " in Table 2) and communications between places were much more costly compared to normal shared memory machines despite the high speed network of UV1000.
Three intervals ns = 10, 100, 1000 were used for load balancing which determined the speed of workload distribution. When the distribution was too slow, the speedup ratio did not scale well (e.g., 3rpr, ǫ = 0.2, with ns = 1000 on the machine (1)). Conversely, small intervals required greater amount of communications and therefore we used ns = 1000 to draw better performance on the cluster where communications were more costly.
There was a large overhead caused by the workers sending a large number of boxes for load balancing when the number of workers was not sufficient against the problem size. Speedups for 3rpr, ǫ = 0.1, using 40 workers or less shows an example of such overheads (Figure 4(d)). Resolving this overhead by suppressing redundant box sends is a part of the future work.

Conclusions
In this paper, we proposed a parallel branch and prune algorithm, based on, non-portfolio, search-space splitting approach. In the experiments, using 256 X10 places (i.e., cores), we achieved speedup factors of as much as 119. We expect that our parallelized solver will be applied to large practical problems, e.g., the robotics problems in [3].