Algorithms for Finding Vulnerabilities and Deploying Additional Sensors in a Region with Obstacles

: Consider a two-dimensional rectangular region guarded by a set of sensors, which may be smart networked surveillance cameras or simpler sensor devices. In order to evaluate the level of security provided by these sensors, it is useful to ﬁnd and evaluate the path with the lowest level of exposure to the sensors. Then, if desired, additional sensors can be placed at strategic locations to increase the level of security provided. General forms of these two problems are presented in this paper. Next, the minimum exposure path is found by ﬁrst using the sensing limits of the sensors to compute an approximate “feasible area” of interest, and then using a grid within this feasible area to search for the minimum exposure path in a systematic manner. Two algorithms are presented for the minimum exposure path problem, and an additional subsequently executed algorithm is proposed for sensor deployment. The proposed algorithms are shown to require signiﬁcantly lower computational complexity than previous methods, with the fastest proposed algorithm requiring O ( n 2.5 ) time, as compared to O ( mn 3 ) for a traditional grid-based search method, where n is the number of sensors, m is the number of obstacles, and certain assumptions are made on the parameter values.

Introduction
A wireless sensor network (WSN) can be used for various tasks such as monitoring the environment and industrial facilities. Each generation of sensors tends to be smaller, more energy efficient, and less expensive than the previous generation. These trends have permitted the development of systems with an extremely large number of sensors. Management of such large networks entails solving problems related to how sensors are deployed, interconnected, and communicate with each other [1].
An important application of sensor networks and WSNs is surveillance [1], which can be conducted using smart networked surveillance cameras (a type of expensive networked sensor device) or simple sensor devices (such as heat, sound, or range sensors). For this type of application, the degree of security provided by the WSN is determined mainly by the positioning of the sensors. A common method for evaluating the adequacy of this type of WSN is to compute the minimum exposure path, which is a path that traverses the guarded region with the lowest level of exposure (Section 3.2) to the sensors. Then, if the resulting level of security is deemed to be insufficient, the existing sensors can be repositioned, additional sensors can be deployed, or both. This paper targets the problem of WSN surveillance in a rectangular region with obstacles and sensors that have sensitivities that vary depending on the distance to the target being detected [2][3][4][5][6][7][8]. The level of security provided by the sensors can be evaluated by finding the minimum exposure path in this target region. This paper proposes a general model and solution for this problem. In the proposed method, the outer limits of the sensing ranges of the sensors and a segment tree data structure are used to compute a feasible area, within which the minimum exposure path must reside. Then, a grid is used within this feasible area to search for the minimum exposure path. Given n sensors, m obstacles, and a few practical assumptions on model parameter values, the computational complexity of the proposed algorithm is O(n 2.5 ), which is a significant improvement over previous approaches, with computational complexity O(mn 3 ).
In addition to computing the minimum exposure path, the proposed algorithm can be used to determine the location (and pointing direction, if appropriate) where an additional sensor can be placed to reduce the vulnerability of the target region. This can easily be done since the process for computing the minimum exposure path includes the determination of the feasible area for this minimum exposure path, and the feasible area is the most vulnerable area [9][10][11]. The main contributions of this paper are as follows. • A general model, which encompasses most of the previously proposed models, is proposed for the minimum exposure problem. • A new algorithm is proposed for the minimum exposure problem. For large n and m values, the proposed algorithm is a significant improvement over previous approaches. • A new near-optimal algorithm is proposed for proper sensor deployment.
The rest of this paper is organized as follows. Section 2 reviews previous related research work. In Section 3, the problem being addressed in this paper is formally modeled and the proposed solution approach is described. In Section 4, the proposed solution and previous possible solutions are compared analytically, and an efficient sensor deployment algorithm is presented. Finally, concluding remarks are provided in Section 5.

Voronoi Diagram and Delaunay Triangulation
The Voronoi diagram [12] is a graph consisting of edges created by drawing dividing lines between each pair of adjacent sensors such that each line is equidistant to the two sensors that it separates, with the line being terminated when it intersects with another line or the outer boundary of the region being monitored. The Delaunay triangulation is the dual graph of a Voronoi diagram. These two graph theory concepts have been found to be useful for computing minimum and maximum exposure paths given certain assumptions on the sensitivities of the sensor used.
Marengoni et al. [1] describe research work on surveillance in a given area referred to as an art gallery. Given that the entire art gallery should be covered by at least one sensor, the authors address the problem of the minimum number of sensors that need to be deployed. The sensors are modeled as disks with infinite ranges. Even though it only considers a two-dimensional (2-D) area, if there are Steiner points, the coverage problem becomes NP-hard [1]. Otherwise, it is NP-complete. To solve the same problem in a 3-D space, the authors of [1] suggest the use of a polynomial approximation solution, which adopts a greedy approach using Delaunay triangulation and graph coloring concepts [12].
Meguerdichian et al. [13,14] suggest an approach that applies a geometrical algorithm to the coverage problem. They assume all n sensors are modeled with identical disks in which the sensitivity decreases as the distance between a point and a sensor increases. They define a maximal breach path as a path in which the distance from the closest sensor and itself is kept at a maximum. They also define a maximal support path as a path in which the distance from the closest sensor and itself is kept at a minimum value. To find the maximal breach path, they construct the Voronoi diagram [12]. After that, they set the weight of an edge in the Voronoi diagram to be the distance between that edge and the nearest sensor. Finally, they find the maximal breach path using a binary search process.
The key value used is the arithmetic average of the biggest and smallest values of the edges. The edges with weights over the key value are inserted in G , and Breadth-First-Search (BFS) is conducted on G . If BFS is successful, the key value is increased following the binary search. Otherwise, the key value is decreased following the binary search. After renewal of the key value, G is reinitialized and edges that are bigger than the key value are inserted in G . This process is continued until the difference between two consecutive key values is less than a tolerance value. The maximal support path is found in a similar manner but with a Delaunay triangulation graph. The maximal breach path can be found by first constructing the Voronoi diagram in O(n log n) time. The number of nodes and edges of the Voronoi diagram is O(n), so each BFS takes O(n) time and the binary search performs BFS in log(n) time. Therefore, the total time complexity becomes O(n log n).
Later, using similar definitions and assumptions, Li et al. [15] suggested computing the shortest path using the sums of the weights of the edges, not the minimum weights of edges as in a maximal breach path. To find the shortest path, they used the Bellman-Ford algorithm. However, if the edge values are nonnegative, Dijkstra's algorithm [16], with computational complexity O(n log n), is faster than the Bellman-Ford algorithm, with computational complexity O(n 2 ).
Hou et al. [17] suggested another method that produces both the maximal breach and maximal support paths using the Gabriel graph [12] and Dijkstra's algorithm. A mathematical proof was introduced to show the correctness of their result.
Mehta et al. [18] defined a breach number tree and introduced Kruskal's algorithm [12] to compute the maximal breach path and maximal support path. Using their method, if the Voronoi diagram with weight edges is constructed first, the paths can be computed in O(|P|), where |P| is the number of edges of the paths, and the values of the paths can be computed in O(1).
In general, given n sensors, the computation of the Voronoi diagram itself can be completed in O(n log n) time because the x and y coordinates of all sensors can be sorted first to determine the pairs of sensors to consider in the creation of the Voronoi diagram. The number of vertices |V| and edges |E| in the Voronoi diagram are O(n), so a shortest-path search algorithm, such as Dijkstra's algorithm, could be used to compute the maximal breach path in time O(|V| 2 ) or O(|E| + |V| log |V|) time using a min-priority queue implemented by a Fibonacci heap [16]. Thus, the overall time complexity for the computation of the maximal breach path using a Voronoi diagram is O(n log n).

Effect of Multiple Sensors Along a Path
Although the Voronoi diagram and Delaunay triangulation are useful concepts, they are found to be insufficient for the computation of the exact minimum exposure paths when the combined effect of multiple sensors are considered. For example, Figure 1 shows that if multiple sensors are clustered together in the lower region of an area while only one sensor is located in the upper region, then the minimum exposure path will be above the set of edges in the Voronoi diagram for this area. A detailed description of this effect is shown in Section 3.2 with an experiment based on the formal definitions of sensors, discover function, and exposure function. For more accurate path computation, the authors of [18][19][20][21] suggested a k-support path problem, and the authors of [19][20][21] solved this problem by implementing a k-th Voronoi diagram that includes the effect of the first to k-th nearest sensors. The authors of [19][20][21] defined the k-th nearest point Voronoi diagram (KNP-Voronoi diagram). Constructing the KNP-Voronoi diagram takes O(k 2 n log n) time and computing the minimum k-support path using Dijkstra's algorithm takes O(k 2 n log n) time. The authors of [22] found the maximal support path in the map, including the effect of not only the n sensors but also m obstacles. The authors of [22] suggested the use of the Constrained and Weighted Voronoi diagram (CW-Voronoi diagram). According to [23], the CW-Voronoi diagram has O(m 2 n 2 + n 4 ) vertices and edges, and constructing the CW-Voronoi diagram takes O(m 2 n 2 + n 4 ) time.
A discover function can be defined as or as a circle with decreasing sensitivity as distance increases, so the discover value is usually defined as discover(x, y) = max s∈S distance(s, (x, y)).
The maximum covered-distance of a path is defined with the definition of a path in Section 3.1, and the path p that maximizes the discover value is called the best coverage path or the best support path. Similarly, the minimum covered-distance can be defined as cover2(p) = min discover(p.x, p.y) and the path p that minimizes the cover value is called the minimum coverage path or the maximal breach path. The authors of [2,3] defined a sensor as a circle and found an optimal path based on intensity and exposure. The authors defined a discover function as the sum of all sensitivity potential function values at (x,y) for all sensors, and an optimal path is the path that minimizes the exposure value. The authors first found an optimal path with one sensor. If the starting point is (1, 0), the goal is (0, 1), there is a sensor s = (0, 0, 1 √ x 2 +y 2 , f s ), and f s is the same as Equation (5) with d = (1, 0) and θ = 2π, then the optimal path is (cos π 2 t, sin π 2 t) and the exposure along this path is π 2 . Based on the optimal path for a single sensor, they suggested a heuristic method for finding an optimal path with several sensors by running the Dijkstra algorithm map divided into horizontal, vertical, and diagonal grid lines. The authors of [4] defined a sensor as a cone and proposed a heuristic optimal path search method with several sensors by running Dijkstra's algorithm in a map divided into horizontal and vertical grid lines.
Later, approaches were proposed for finding an optimal exposure path. In [24,25], for a single sensor map, an optimal exposure path was calculated by solving the Euler-Lagrangian equation in polar coordinates, while setting ρ and φ as independent variables. The result is exactly the same as [2,3]. The authors of [26] used a similar method, but the optimal exposure path was calculated by solving the Euler-Lagrangian equation in Euclidian coordinates with x and y used as independent variables. For a map with multiple sensors, the authors of [24][25][26] used timing segmentation. They calculated the optimal direction for a specific position for an object, moved it slightly, calculated it again, moved it slightly, and so on until the object arrived at the goal. The authors of [27] set y as a function of x and solved the Euler-Lagrangian equation with x used as an independent variable. For multiple sensors, the authors of [27] considered only neighboring sensors and found an approximated path, which was a localized minimum exposure path.
The sensitivity potential function can be modeled as varying according to the target velocity, d dt p(t) = ( dx dt , dy dt ) [24][25][26]. In addition, the probability of detection will decrease if the target accelerates. Consider a sensor modeled as s = (0, 0, where f s is the sensitivity of the sensor, which varies depending on the sensing location and the velocity of the target. Then the exposure can be modeled as follows.
and the Euler-Lagrangian equation becomes and the solution is which is not feasible. Moreover, the Euler-Lagrangian equation does not provide a solution when the discover function does not have a continuous first partial derivative. Creating a grid on the map and running the Dijkstra algorithm with that grid is an alternative practical solution for general maps. The authors of [5] divided a map with square grids and introduced Gaussian noise in grid edges. After the values of grid edges were set, they found an optimal path with the Dijkstra algorithm. The authors of [28] defined a mathematical model of physarum, used physical features of physarum and found a path using these features. The authors of [6][7][8] not only set the value of grid edges as probabilistic values but also defined an exposure value in a probabilistic manner as follows.
discover : Furthermore, the multiplications can be changed to additions by converting the values to their logarithms. Then, additions can be used to calculate the grid edge values in the integral.
After that, they used the Dijkstra algorithm to find the shortest path. The authors of [29] use both the Voronoi diagram and the exposure concept. They used the Voronoi diagram but then corrected the path by considering the exposure to neighboring sensors. The authors of [30] then proposed an improvement by defining intensity as the sensing value of the nearest sensor and then integrating the effect of all sensors.

Sensor Deployment
The authors of [2,3,6,8,13,14,17,[30][31][32] defined sensors that have circle-shaped sensitivity functions. The authors of [13,17] deployed sensors in the middle of the least weighted edge of the Voronoi diagram. The authors of [2,3] conducted an experiment that compared random deployment and deterministic deployment. In deterministic deployment, they deployed sensors in the middle sections of their grid edges, which can be formed into triangle or hexagon shapes by changing the definitions of sensors and the discover value. Reference [13] also divided the target map into grids and formulated the one-level coverage problem, and then proposed a heuristic algorithm called Tabu Search (TS). The authors of [30] deployed sensor networks based on Average Linear Uncovered Length (ALUL) and Minimum Linear Uncovered Length (MLUL) definitions. They divided the map of grid edges and selected the position of the additional sensor that set the ALUL or MLUL to be the maximum. For MLUL, they used the Voronoi diagram. The authors of [6,8] deployed some of the sensors randomly and checked whether the minimum exposure satisfied the desired level. If not, more sensors were added and the previous step was repeated. This was continued until the minimum exposure satisfied the desired level. Reference [14] deployed an additional sensor in the middle of the grid edge square in which the summation of the edges was the smallest. Reference [31] deployed sensors randomly the first time and then let each sensor have a virtual repulsive force to other sensors so that they pushed each other apart. The authors of [33] also deployed sensors with the virtual force for minimizing the required number of mobile sensors in sweep coverage under the energy constraint. The authors of [32] presented an algorithm for finding locally-optimal coverage in the minimum unit-disk cover problem through a Lagrangian minimum approach with O(ln 4 ), where l is the number of iterations and n is the number of locations that must be included.
The authors of [9][10][11][34][35][36][37] defined sensors that have cone shaped sensitivity functions, which is a more general form than a circle. The authors of [34][35][36][37] deployed sensors in a region with obstacles. The authors defined a mathematical model for sensors in which the sensitivity extends to the boundary of the region or obstacles. The authors of [34,35] deployed sensors with a greedy algorithm. Reference [36] divided a region into grids and suggested a solution based on solving the integer nonlinear programming problem. As this becomes more difficult when the problem size gets larger, the authors of [37] suggested an immune-based algorithm. The authors of [38] suggested a barrier sweep coverage problem for covering not a whole area but finite critical curves with sector type sensors. They proved the barrier sweep coverage problem is an NP-hard problem and proposed an approximate algorithm.
The authors of [9,10] defined the most general form of the sensor, which has its own sensing radius and directional vector. They used virtual repulsive force and calculated the directional vector by virtual force. The authors of [11,39] applied Particle Swarm Optimization (PSO) to the coverage problem with this sensor model, and Reference [40] suggested Hybrid Genetic Particle Swarm Optimization (H-GPSO), which is a mixing model of the models in Reference [28] and [11]. The region considered in [9][10][11] does not include obstacles.

System Model
A general model is presented for the computation of the minimum exposure path problem. This model will attempt to encompass almost all of the previous models used while addressing their shortcomings.
We assume that the region to be guarded is a 2-D rectangle of dimensions B X and B Y . n sensors s i ∈ S are positioned within this region to detect trespassers, but the approach used here can easily be extended to a 3-D arbitrarily-shaped guarded region by replacing 2-D vectors to 3-D vectors. Each sensor s i ∈ S is positioned at s i,pos = (x i , y i ) and has a sensitivity potential function f s i ( x pos ) that denotes the strength of the signal detected at sensor s i due to a trespasser located at x pos = (x, y). Depending on the function f s i ( x pos ), we can model a sensor for which an event within a radius of (x i , y i ) is always detected whereas all events farther than the radius from s i remain undetected (the commonly-used solid disk model). More realistically, an f s i ( x pos ) function that decreases with the distance to x pos = (x, y) can be used to model a sensor in which the sensitivity to an event diminishes as the distance to the event is increased. Finally, this type of general definition for f s i ( x pos ) can also be used to model a sensor that can monitor events within a cone shape, such as would be the case with a camera with a limited field of view.
The presence of a set W of impassable obstacles is assumed, such as pillars, facilities in a factory, or a great pile of luggage. For an obstacle w ∈ W, a function f w exists that describes the obstacle as w = {(x, y) (x, y) ∈ X × Y and f w ( x pos ) > 0}. The function Boundary(area) is used to denote the boundary of a given area. For example, given f w as below, obstacle w is a rectangular area (Section 3.3). Thus, If the sensing area occupies the back of the obstacle, it means that the sensor has a field of view on the back of the obstacle, so it implies the situation should be a 3-D model rather than a 2-D model. Therefore, it is assumed that there is no field of view behind the obstacle from the 2-D sensor.
For example, Figure 2 shows a university dormitory area with four CCD (Charge Coupled Device) cameras and one dormitory building. A thief is at the start location, and the thief wants to get to the goal location while avoiding the sensing area of CCD cameras as much as possible. As the dormitory building is higher than any camera height, the camera's field of view cannot reach behind the building (Figure 2). The sensitivity potential function for a sensor s ∈ S is assumed to have the following general form, where η ∈ R is the noise potential, d s is a directional vector, and θ s 2 is the angle of the fieldof-view; K s , l 1 , and l 2 are sensor-dependent constants (Figure 3). The given coordinate x pos and position of sensor s i s i,pos are determined according to the dimension of the guaranteed region, which is (x, y) and (x i , y i ) for 2-D ( Figure 3a) and (x, y, z) and (x i , y i , z i ) for 3-D (Figure 3b).  In addition, we adopt the definitions for coverage and sensing area used in [41] and generalized the definitions with obstacles. To prevent a target from penetrating into obstacles, we define the sensitivity potential in obstacles to be positive infinity. In X × Y, we should construct the appropriate data structure using W and S. The appropriate data structure is one for which, when a coordinate is given, the information about which sensor functions or obstacles are covered is returned. We assume that sensors cannot see the target when it is obscured by obstacles. As an example, this assumption and Definition 1 results in sensing areas S1 to S11 when the sensors and an obstacle are arranged as in Section 3.3, which is geometrically simplified from the realistic scenario in Figure 2.
Finally, the area can be defined formally. The position where the target starts to move, or the entrance of the map, is denoted as (x in , y in ) ∈ X × Y, and the position where the target stops moving, or the exit point of the map, is denoted as (x out , y out ) ∈ X × Y. The map M is denoted formally as (X, Y, S, W, (x in , y in ), (x out , y out ), η), where η ∈ R is the noise potential.

Discover and Exposure Functions
The target wants to move from (x in , y in ) to (x out , y out ) with the least chance of being discovered by a sensor. To formally find the movement path of the target, we can consider a real analytic function p from the time interval to a specific point in the map, p : R → X × Y. However, the paths we consider start at (x in , y in ) and end at (x out , y out ). To avoid the infinite velocity value that results from the Euler-Lagrangian Equation (2), we set dp(t) dt = 1. Finally, the set of paths is formally denoted as P = {p , p is real analytic, the initial point of p is p start , and the final point of p is p goal }. A path p is an element of P.
The calculations using Equations (1) and (3) to set the optimal path are similar. Eventually they calculate the discover value by summing several functions and then calculate exposure as the integral of the discover function; the minimum paths are the same because the exponential function is increasing. For this reason we select Equation (1) for our target calculations.

Definition 2.
A discover value about specific position x pos is discover( x pos ) = ∑ s∈S f s ( x pos ). An exposure value about specific path p is exposure(p) = We also set the path that has the smallest exposure value of map M as p min (M). The exposure(p min (M)) is the degree of vulnerability of a region, which can be defined mathematically as shown below.
Although the Voronoi diagram and Delaunay triangulation are useful concepts, they are found to be insufficient for the computation of the exact minimal exposure paths when the combined effect of multiple sensors is considered. For example, if multiple sensors are clustered together in the corner of the map with weighted values (Figure 4), then the minimal exposure path from the Voronoi diagram is not actually the minimum. To verify that, we set an experiment. From the definition of the sensor (Equation (5)), we set θ s = 2π, K s = 10,000, l = 2, and set the noise η = 1. If the sensors are deployed following Figure 4, the resulted minimum exposure path from the Voronoi diagram usually represents more exposure value even than the 50 × 50 grids diagram ( Figure 5).

Feasible Paths and Areas
The proposed method reduces the time required to find the minimum exposure path. To accomplish this, we need to define additional sets and functions. Section 4 presents the actual algorithms and proofs that use these sets.
p min (N, A) is the path found by running the Dijkstra algorithm following N grid edges in area A.p min (N, A) approaches p min (M) as N goes to positive infinity, However, as N increases,p min (N, A) becomes more accurate but the number of computations increase. To reduce the number of computations with the same level of accuracy, we decided to limit the area of interest in the target map, and this area is called the feasible area,Ã. To findÃ, the skeleton of the feasible area, feasible paths, should be defined. A set of paths following the boundaries of sensing areas is denoted asP = {p(t) p ∈ P and p ∈ Boundary(whole sensing areas)}, and p min (M) = arg min p∈P exposure(p).

Minimum Exposure Path
Definitions for discover and exposure functions were presented in Section 2.2. References [2][3][4][5][6][7][8] tried to find p min (M) by partitioning the map area using edges and weighting each edge using the exposure value along the edge. Finally, they ran the Dijkstra algorithm while following the grid edges (Algorithm 1).
This paper proposes a new, more efficient method for finding the minimum exposure path, as shown in psuedocode format in Section 4.1. In the proposed method, we want to reduce the time consumption by reducing the time required to calculate each grid edge, and reducing the number of grid edges to be considered. To reduce the time required to calculate each grid edge, we added the step "Construct sensing area data structure" and set the grid edge values "based on an efficient data structure" (Algorithm 2). For additional time reduction, we added the step "MakeÃ" and "Run the Dijkstra algorithm with grid edges inÃ" (Algorithm 3). To "Run the Dijkstra algorithm withP" (Figure 8), the vertices and edges can be like those shown in Figure 8. The Dijkstra Algorithm can be used to find the shortest path using the vertices and edges shown. We will show that Algorithms 1 and 3 return the same result, but that Algorithm 3 is faster than Algorithm 1. Proof of Theorem 1. The sensing area data structure uses two concepts: a doubly-connected edge list and a segment tree. The doubly-connected edge list can be used to determine which sensors cover a specific region in O(d 1 ) time if a boundary edge is given, where d 1 is the maximum number of sensors in a sensing area. The segment tree can be used to find a boundary edge in O(log(m + n)) time if a point is given. When taken together, the list of sensors that cover a region can be found in O(log(m + n) + d 1 ) if a point is given. This is much faster than O(m + n).
When a sensor is included in the map, intersections and the leftmost point and the rightmost point become vertices. A single line segment bifurcates into two edges if they have different exposure values. This rule has two exceptions: when the line segment is outside the map boundary and when it is inside an obstacle, because a target cannot go outside the map boundary or inside an obstacle (Figure 9). The edges created in this manner are stored in a segment tree ( Figure 10). For each sensing area, information about edges, vertices, and the set of sensors that affect the sensitivity in that sensing area are stored.
The number of vertices, edges, and areas is O((m + n) 2 ). Therefore, the sensing data structure can be constructed in O((m + n) 2 log(m + n) 2 ) = O((m + n) 2 log(m + n)) time.  In Algorithm 3,Ã is the gray region in Figure 7. Grids are calculated only in the gray region in Figure 7. Let N X be the number of lines vertically, and N Y be the number of lines horizontally. The grid lines make a grid square. Then ) and the time complexity to calculate all edge values is O(N 2 X (d 2 + log(m + n) + d 1 )) and the time complexity to find the optimal path is O(N 2 X log N 2 X ) = O(N 2 X log N X ). Therefore, the total time complexity is = O((m + n) 2 log(m + n) + N 2 X (d 2 + log(m + n) + log N X + d 1 ). Furthermore, by applying 0 ≤ k ≤ d 1 ≤ n and 1 ≤ d 2 ≤ (m + n), the total time complexity of Algorithm 2 becomes O((m + n) 2 log(m + n) + N 2 X (m + n + log N X )) Now we show why the "Construct sensing area data structure" step should be included. Previous work [2,3,[24][25][26][27] did not include this step. If this step is ignored, all sensors and obstacles must be considered for all grid edges. The time complexity to cal-culate values for all edge values is O(N 2 X mn), and to findp min (N, . Therefore, the total time complexity of Algorithm 1 will be O(N 2 X mn + N 2 X log N X ) = O(N 2 X (mn + log N X )). This is obviously bigger than O(N 2 X (m + n + log N X )). To find a more accurate path than the path found with a Voronoi diagram, N X ≥ O(m + n). This makes O(N 2 X (mn + log N X )) > O((m + n) 2 log(m + n)). Therefore, including the "Construct sensing area data structure" step reduces the necessary time complexity. Proof of Theorem 3. At first, the path can be filtered if the sensitivity function of each sensor is convex downward. To simplify the proof, we consider a fractional sensitivity potential model modeled by Equation (5) (Figure 11). For an arbitrary convex downward function y = g(x), g (x) ≥ 0, and g (x + L) ≥ g (x) for all x and L > 0. We set h(x, L) = g(x + L) − (g(x) + Lg (x)), from which it follows that h(x,L) ∂x = g (x + L) − g (x) + Lg (x) ≥ 0 and h(x, 0) = 0. Therefore, h(x, L) ≥ 0 for all x and L > 0.
Next, in a map, if we apply For this reason, any sensitivity potential projected in any direction is convex downward.
No matter how many convex downward functions are added, the added function is also convex downward. For these two reasons, outside ofÃ, no path has an exposure value lower thanp min (N,Ã).Ã always containsp min (N, M). Moreover,Ã contains p min (M) because the maximum distance betweenp min (N, M) and p min (M) is less than the size of a side of a cell in the grid.  The number of areas is O((m + n) 2 ). Accordingly, the number of vertices and edges for "Run the Dijkstra algorithm with Boundary(sensing_areas) and configureÃ" are |V| = O((m + n) 2 ), |E| = O((m + n) 2 ). Then the Dijkstra algorithm takes O((m + n) 2 log(m + n) 2 ) = O((m + n) 2 log(m + n)) time. Each of the result paths is stored including the sensor information for the path. In step "DivideÃ with grid edges and set the values from sensing area data structure", grids can be configured through the sensor information. Therefore, those additional parts take O((m + n) 2 log(m + n)) time. This is the same as the step "Construct the sensing area data structure", so this step does not increase the time complexity. For those reasons, the additional computations used in Algorithm 3 do not increase the time complexity.
Although the above analysis has shown the time complexities of Algorithm 1 through Algorithm 3 in as exact a manner as possible, these time complexity equations can be expressed in simpler forms by adopting several practical assumptions.
First, for the reasonable coverage of the map, the number of sensors n should be more than the number of obstacles m, and sensors almost cover the whole map. Then, the time complexity of Algorithm 2 becomes O(n 2 log n + N 2 X (n + log N X )). Next, the shorter given path is, the smaller exposure of the path is. If a feasible path is about a straight line from the start in the lower left corner to the goal in the upper right, the area ratio ofÃ to whole map is about , where S av is the average sensing area of sensors ( Figure 7). According to the first assumption, nS av = N X B X N Y B Y . This induces the number of vertices and edges in "DivideÃ with grid edges and set the values from sensing area data structure" as Then, the time complexity to calculate all edge values is O( n) and the time complexity to find the optimal path is O( n log N X ). Therefore, the time complexity to find the optimal path is O(n 2 log n + N 2 n log N X ). Finally, by considering the computation limits, the number of vertical grid lines N X can be set to O(n), which derives the time complexity of Algorithms 1-3 as, respectively, O(mn 3 ), O(n 3 ), and O(n 2.5 ).

Deployment of Additional Sensors
Motivated byÃ, to handle the vulnerable area, additional sensors should be deployed in a manner that increases exposure(p min (M)) as much as possible. The candidate points are near the locally maximum curvature points due to the definition of the sensitivity potential function. The points at which the curvature are locally maximum are stored in the setL. An example of this set of points is shown asL = {(x 1 ,ỹ 1 ), (x 2 ,ỹ 2 ), (x 3 ,ỹ 3 ), (x 4 ,ỹ 4 ), (x 5 ,ỹ 5 )} in Figure 12.
For a point at which the curvature is a local maximum, such as (x j ,ỹ j ), where 0 ≤ j ≤ |L|, the normal vector − → N j can be computed. − → d j is set to the opposite direction of − → N j , so − → d j = − − → N j because the directional vector and normal vector are unit vectors. Let (x j ,ȳ j ) be the position of s for the outer curve of Boundary(Ã), with the sensing area of s being tangent to each other ( Figure 13).
To compare the points inL, the maximum increment coverage value should be selected. However, the increment value is only considered in the line between (x j ,ȳ j ) and (x j ,ỹ j ) because the error is not considerable. Let L j be the distance between (x j ,ȳ j ) and (x j ,ỹ j ) ( Figure 14). By following variable ξ j between [0, L j ], the position (x j , y j ) = (x j ,ȳ j ) + ξ j L j − → d j varies. The outer interval of [0, L j ] is not considered because the area of (Ã ∩ sensing area of s) by setting (x s , y s ) = (x j , y j ) becomes smaller in the outer of the interval of the position of s. Therefore, the approximate maximum increment coverage value for (x j ,ỹ j ) is denoted as incr j , and the corresponding ξ j is denoted as ξ j,max . incr j = max 0≤ξ j ≤L j Ã ∩sensing area of s, where (xs,ys )=(x j ,y j ) f s (x, y)dxdy, ξ j,max = arg max 0≤ξ j ≤L j Ã ∩sensing area of s, where (xs,ys )=(x j ,y j ) f s (x, y)dxdy.
For each (x j ,ỹ j ), incr j is calculated and all locally maximum curvature points are compared based on the incr j value. The index for which incr j is the maximum value is denoted as j max . j max = arg max 0≤j≤|L| incr j .
From j max , the undefined components of the additional sensor s are now defined, (x s , y s ) = (x j max , y j max ), and − → d s = − → d j max . After the sensor s is included in M, exposure(p min (M)) will be increased deterministically because the approximate position for the maximum increase is considered (Algorithm 4). Calculate (x j ,ȳ j ), and L j 7: Calculate incr j and ξ j,max 8: end for 9: Calculate j max , and fill the undefined components of s 10: Update sensing area with s 11: return M

Conclusions
This paper has addressed the general problem of WSN surveillance security by proposing a minimum exposure path search problem and a sensor deployment problem. Sensors, obstacles, and the target surveillance area have been modeled, and computationally efficient solutions have been proposed. For the minimum exposure path problem, the sensing limits of the sensors and the obstacle locations are used to form a "feasible area", within which the search is constrained. Vertical and horizontal grid lines are then used in this feasible area in order to search for the minimum exposure path in a computationally efficient manner. The proposed algorithm has computational complexity O(n 2.5 ), a significant improvement over previous approaches, with complexity O(mn 3 ), where m is the number of obstacles and n is the number of sensors.


Local Search with a SAT Oracle for Combinatorial Optimization

NP-hard combinatorial optimization problems are pivotal in science and business. There exists a variety of approaches for solving such problems, but for problems with complex constraints and objective functions, local search algorithms scale the best. Such algorithms usually assume that finding a non-optimal solution with no other requirements is easy. However, what if it is NP-hard? In such case, a SAT solver can be used for finding the initial solution, but how can one continue solving the optimization problem? We offer a generic methodology, called Local Search with SAT Oracle (LSSO), to solve such problems. LSSO facilitates implementation of advanced local search methods, such as variable neighbourhood search, hill climbing and iterated local search, while using a SAT solver as an oracle. We have successfully applied our approach to solve a critical industrial problem of cell placement and productized our solution at Intel.

Loosely speaking, an instance of a combinatorial optimization problem deals with the minimization of an objective function over a finite set, subject to feasibility constraints (or, simply, constraints). The set of all elements satisfying the constraints is referred to as the set of feasible solutions (or, simply, solutions). In this paper, we focus on solving any problem, which can be expressed as a constraint optimization program (COP) [2]. Arguably, the vast majority of combinatorial problems, encountered in practice, fall under this category.
Many important combinatorial problems are NP-hard. For such problems, various algorithmic strategies have been devised, including complete methods, such as branchand-bound and dynamic programming, and incomplete methods, such as greedy algorithms and local search. Each such method imposes requirements on the mathematical properties of the problem with a consequent limit on the scope of applicability.
Local search algorithms stand out from the rest in that they impose relatively mild constraints on the type of the problem to be addressed, thus providing a wide scope of applicability. Furthermore, they seem to scale better with input size relative to complete algorithms [24]. This makes local search algorithms an attractive choice. However, local search algorithms may return a low-quality solution or no solution at all, given a problem for which the mere task of finding a feasible solution is NP-hard. Henceforth, we shall refer to such problems as NP-Hard-Feasible problems. This paper introduces the Local Search with SAT Oracle (LSSO) methodology, that is, local search algorithms which use a SAT solver (or a SAT-based optimization algorithm; details appear later) as an oracle. A key advantage of our approach is that it can handle problems with complex constraints and objective functions. In particular, it can handle NP-Hard-Feasible problems.
To see how SAT solvers might be useful, consider the basic version of a local search for an optimal solution. At the beginning, the local search generates an initial solution and sets it as the current solution. Then, it enters a loop. In each iteration, it looks for a solution with a lower value of the objective function within a neighbourhood of the current one. If such a solution is found, it is set to be the current solution, and the execution resumes. Otherwise, the algorithm terminates and returns the current solution.
A key component of local search algorithms is the neighbourhood function, which assigns to each feasible solution a subset of feasible solutions, called its neighbourhood. Ordinarily, a neighbourhood of the current feasible solution comprises a set of solutions which can be obtained from the current solution by applying a small collection of feasibility-preserving perturbations to its combinatorial structure. A key concern is ensuring that neighbourhoods: (i) are polynomially searchable, and (ii) contain high-quality solutions. However, meeting both requirements might be challenging, since polynomial searchability implies that neighbourhoods should be small, and hence less likely to contain high-quality solutions. In addition, in the case of NP-Hard-Feasible problems, it is not clear how to achieve polynomial searchability, since a search should, in particular, be able to find a feasible solution, which is NP-hard.
Our main idea is to let the SAT solver both find an initial solution and conduct the neighbourhood search. The designer can now define feasibility constraints and neighbourhoods declaratively, that is, by a set of SAT constraints. The designer has more freedom to choose neighbourhoods, which need neither be small, nor contain only solutions close to the current solution. This is because the search of the now complex and possibly large neighbourhoods is entrusted to SAT solvers, constructed precisely to efficiently search large complex subspaces. Our approach lends itself to implementations of advanced local search variants, such as variable neighbourhood search, hill climbing and iterated local search [29].
An important feature of our algorithms is that they are anytime. Recall that an anytime algorithm is expected to return a valid solution even if interrupted. An anytime algorithm for an optimization problem is expected to find an improving set of solutions. The anytime property is essential for industrial application, since it allows the user to get an approximate solution even for very difficult instances [14,15].
We demonstrate the usefulness of our approach by solving hard industrial instances of the NP-Hard-Feasible cell placement problem. Cell placement is one of the most important problems in VLSI automation [28]. Its most basic version concerns placing without overlap a set of rectangles on a grid, while minimizing the occupied area. In reality, the problem is more complex. Our approach has been successfully productized at Intel.
The rest of this paper is organized as follows: Sect. 2 provides the necessary background. Sect. 3 introduces our LSSO methodology. Sect. 4 shows how to solve placement with LSSO. Sect. 5 presents the experimental results. Sect. 6 concludes our paper.

Background
This section provides some background. Sect. 2.1 is an overview of COP. Sect. 2.2 describes the cell placement problem and shows how to reduce it to COP. Sect. 2.3 discusses how one can solve a COP using a SAT-based bit-vector solver. Sect. 2.4 reviews local search.

Constraint Optimization Program (COP)
This work presents a new methodology for solving a wide class of combinatorial optimization problems, which can be expressed as a Constraint Optimization Program, shown in Def. 1.
Definition 1 (Constraint Optimization Program (COP) [2]). A constraint optimization program is a tuple (X , D, C, Ψ ) where: . . x n } is a finite set of variables often referred to as decision variables. 2. D = {D 1 . . . D n } is a corresponding set of finite domains. Without loss of generality, each D i is assumed to be a closed bounded interval of non-negative integers.
is an objective function to be minimized.

The Cell Placement Problem
Cell Placement (Placement) is a major stage in the VLSI design cycle [8,16]. The input of the cell placement problem comprises the following components: 1. A rectangular grid region of M rows and N columns, on which the cells are to be placed. Row/column line numbering starts at 0 and ends at M /N , respectively. 2. A finite set C of rectangular cells. The width and the height of each cell c ∈ C are assumed to be positive integers, denoted by c width : 0 < c width ≤ N and c height : 0 < c height ≤ M , respectively. 3. A set R of forbidden rectangular regions. A forbidden region r ∈ R is specified by 4 numbers r west , r south , r east and r north (where, 0 ≤ r west , r east < N ; 0 ≤ r south , r north < M ; r east > r west ; r north > r south ), denoting the leftmost column line, bottom row line, rightmost column line, and top row line, respectively. 4. A finite set I of nets, each consisting of a non-empty subset of cells. The nets may (and usually do) intersect.
We are interested in feasible placements, that is, placements in which no cell overlaps other cells or forbidden regions. Given a feasible placement, we define the size of a net n ∈ I as the perimeter of the box bounding its placed cells. We define the size of the placement as the sum of the sizes of the nets. We are required to find a feasible placement of a minimal size. An example is shown in Fig. 1.
In industrial practice, there may be additional industrial requirements, such as aligning some of the cells, enforcing parity constraints (i.e., the user might require the y coordinates of some of the cells to be either even or odd) [19], ensuring a minimal distance between some of the cells and others. We omit further details due to IP considerations.
Placement is NP-Hard-Feasible, since the NP-complete bin packing problem can be reduced to the decision version of the placement problem [10].
2.2.1 Constraint Optimization Program for Cell Placement. We show how to construct a COP for the cell placement problem. For each cell c ∈ C, let c west and c east denote its leftmost and rightmost column respectively, and c south and c north denote its bottom and top row. Strictly speaking, it suffices to use c west and c south as the COP's independent variables, but it is convenient to use c east and c north as syntactic sugar for c west + c width and c south + c height , respectively. The COP looks as follows: 3. Feasibility constraints: (a) Each cell c is placed wholly within the grid region: For every pair of cells c i , c j , such that i < j, there is no overlap: For every pair r, c of a forbidden region r and a cell c, there is no overlap: (d) Constraints representing any additional industrial requirements. 4. Objective function Ψ : for every net n ∈ I, let n denote its size. We have:  [16]. A solution is shown for the problem of placing five cells c1, c2, c3, c4 and c5 of sizes 4×1, 4×3, 2×2, 2×4 and 1×5 respectively, on a grid with M = N = 8. There are three nets: n1 = {c1, c3, c5}, n2 = {c2, c3} and n3 = {c2, c4} (without any forbidden regions). The bounding boxes of the nets are B1, B2 and B3, respectively. The sizes of the nets, comprising the perimeters of the bounding boxes, are 20, 18 and 20, respectively. The overall placement size is 20 + 18 + 20 = 58. The solution is an optimal one.

Solving COP with SAT
A COP can be solved with various types of solvers [2]. In particular, it is possible to solve a COP by reduction to a series of SAT solver invocations through bit-vector reasoning as explained below.
2.3.1 Bit-vector Solving and SAT. We start with reviewing the basic terminology, related to SAT solving. A literal l is a Boolean variable v or its negation ¬v. A clause is a disjunction of literals. A formula F is in Conjunctive Normal Form (CNF) if it is a conjunction (set) of clauses.
A SAT solver [4] receives a CNF formula F and returns a satisfying assignment (aka, model or solution), if one exists. In incremental SAT solving under assumptions [5,18], the user may invoke the SAT solver multiple times, each time with a different set of assumption literals (called, simply, the assumptions) and, possibly, additional clauses. The solver then checks the satisfiability of all the clauses provided so far, while enforcing the values of the current assumptions.
A bit-vector variable (bit-vector) of width n = |B|, B = {v n , v n−1 , . . . , v 1 }, is a sequence of n Boolean variables, called bits. Bit v 1 is the Least Significant Bit (LSB) and v n is the Most Significant Bit (MSB). A Boolean constant is either ⊥ (0) or (1). A bit-vector constant is a bit-vector (BV), each one of whose bits is substituted by a Boolean constant. A bit-vector term is either a bit-vector, a BV constant, or a result of applying an operator which returns a bit-vector (for example, BV addition, if-thenelse, concatenation) over other terms and atoms. An atom is either a Boolean variable, a Boolean constant or a result of applying an operator, which returns a Boolean (for example, = or unsigned-less-than), over BV terms and atoms. A bit-vector formula (also known as a bit-vector constraint) is recursively defined to be either an atom, a negation of a bit-vector formula, or the result of applying the Boolean operator ∧ or the Boolean operator ∨ over two or more bit-vector formulas. See [3,12] for a rigorous description of the BV language. A BV solver decides the satisfiability of BV formulas.
A BV formula F is satisfiable iff it has a model, that is, an assignment of BV and Boolean constants to their corresponding BV and Boolean variables, which satisfies F . In this paper, BV constants are interpreted as unsigned numbers, and BV comparison operators are interpreted as unsigned. For example, given a bit-vector B = {v 3 , v 2 , v 1 }, the formula F = B < 2 has two models µ 1 : µ 1 (B) = 0 and µ 2 : µ 2 (B) = 1.
All the algorithms presented in this work are assumed to use the so-called eager BV solver [6] which, following some preprocessing, translates the input BV formula to an equisatisfiable formula in CNF and solves it with a SAT solver. Thus, we will use the notions of BV solving and SAT solving interchangeably. We also assume the BV solver to have the same incremental API as a SAT solver.
Since the variables in a COP have finite domains, both the variables and the constraints of a COP can be easily expressed as BV variables and BV constraints.
In particular, in the COP constructed for the cell placement problem in Sect. 2.2.1, the variables and the constraints can be expressed as BV variables and constraints as follows: For each cell c, we define four bit-vectors: c west and c east of width log N as well as c south and c north of width log M . All the constraints in our COP involve these bit-vectors and can be expressed in terms of operators and relations available in the BV language [3]. Specifically, we implement min and max operators using a series of if-then-else operators. In addition, for every operator, we zero-extend the widths of the operands and the resulting bit-vector to prevent an overflow, whenever required.
Reducing the constraints of a COP to a BV formula and invoking BV solver suffices to find one non-optimal solution. However, for solving the optimization problem by reduction to BV, one needs an extension of BV solving to optimization. 1

Extending
Bit-vector Solving to Optimization. One can extend bit-vector solving to the so-called Bit-Vector Optimization (OBV) [19] as follows: A model µ of a BV formula F is T -minimal, for a given bit-vector T , iff µ(T ) ≤ ν(T ) (where the comparison is unsigned) for every model ν of F . Given a BV formula F and a term T = {t n , t n−1 , . . . , t 1 } in F , where T is called the optimization target (or, simply, the target), Bit-Vector Optimization (OBV) is the problem of finding a Tminimal model of F . The bits of the target T are referred to as the target bits.
Translating our placement COP to OBV is straightforward. We have already shown how to translate the constraints. The optimization target is constructed in the same way as the objective function Ψ is constructed in the COP.
How can one solve OBV in practice? First, one can use the following simple anytime Linear Search algorithm, implemented on top of an incremental BV solver [16,27]: 1: solver.Assert(F ); µ := solver.Sat() assert F and find the first solution 2: while µ is a solution do while there is still a solution 3: solver.Assert(T < µ(T )) block all the solutions with cost ≥ µ(T ) 4: µ := solver.Sat() can we improve our solution?
5: return µ µ is guaranteed to be T -minimal Another anytime algorithm to solve OBV is the following binary search-based algorithm, called OBV-BS [9,19]: 1: solver.Assert(F ); µ := solver.Sat() assert F and find the first solution 2: i := n i is the current bit number, initialized to the MSB 3: while i ≥ 1 and µ(t i ) = ⊥ do fix to ⊥ the MSBs, assigned to ⊥ in µ 4: solver.Assert(¬t i )

5:
i := i − 1 after the loop, i will point to the first target bit, assigned 6: while i ≥ 1 do Check one-by-one, if we can flip the remaining target bits to ⊥ solver.Assert(t i ); i := i − 1 t i cannot be flipped to ⊥, so we fix it to 12: return µ We have successfully applied OBV-BS for solving the problem of fixing an existing placement [19], closely related to the generic placement problem, we are exploring in this work. However, both Linear Search and OBV-BS failed to scale to industrial instances of our current problem of finding an optimal placement from scratch (with Linear Search scaling somewhat better than OBV-BS).
Recently, we have introduced the so-called Polosat anytime algorithm [16], which can be used instead of the standard SAT solver inside Linear Search (and other SAT-based anytime optimization algorithms) to make it substantially more scalable. The idea behind Polosat, shown below, is to simulate local search using a SAT solver. We use the strictly-monotone version of Polosat [16], which assumes the availability of the so-called Boolean observable variables (observables) Obs, that is, a set of Boolean variables on which the objective function depends (for placement, the observables might comprise the bits of the bit-vectors, representing the sizes of the nets, for every net). Polosat is carried out by getting a model µ and then trying to improve it by repeatedly flipping observables, which have not been assigned ⊥ in previous models: 1: function SOLVER.POLOSAT(assumptions) Require: Target bit-vector T is available; Observables Obs are available. return µ To combine Polosat into Linear Search, it is sufficient to replace solver.Sat invocations by solver.Polosat invocations in the code. 2 We have shown in [16] that replacing plain SAT invocations by Polosat invocations in Linear Search makes our placement tool substantially more scalable. We reaffirm this result in Sect. 5.
Yet, despite the significant progress we had witnessed when applying Polosat, we found that combining Polosat into Linear Search is still insufficient for solving a variety of complex real-world instances of our industrial placement problem. This empirical challenge lead us to develop our LSSO methodology, presented in this paper. As we shall see, combining LSSO and Polosat makes our tool considerably more scalable, while the methodology itself is generic and can be applied to solving a wide range of optimization problems.

Local Search Algorithms
Local search strategies [1] are a collection of algorithmic templates. An algorithmic template specifies the main flow of an algorithm, but leaves some details unimple-mented. By implementing these details for a specific problem, one obtains an algorithmic solution for that problem.

Basic Local Search Strategy.
The basic strategy generates an initial feasible solution and sets it as the current solution. Then, it enters a loop. In each iteration, it looks within a neighbourhood of the current solution for a feasible solution with a lower value of the objective function. If one is found, it is set to be the current solution. Otherwise, the algorithm is terminated returning the current solution. Note that this version is guaranteed to stop; it does so, when it reaches a local minimum of the objective function with respect to the neighbourhood used.
To turn this algorithmic template into a complete algorithm, one has to implement the following problem-dependent items: (i) A procedure for generating an initial feasible element. (ii) A neighbourhood function assigning to each solution a subset of solutions. (iii) An algorithm for searching the neighbourhood for a better solution.

Neighbourhood Functions.
A key factor, which affects both the complexity of the search and the quality of the resulting solution, is the selection of a neighbourhood function. In theory, the selection ought to depend on a mathematical analysis of the structure of the feasible set and the profile of the objective function. For complex problems, such an analysis is usually beyond reach. The classical approach to neighbourhood definition is based on the following problem-independent general principles: 1. Drawing on an analogy to optimization algorithms in the continuous case (such as gradient descent or line search), a neighbourhood should be so defined as to make its elements "close" to the current solution. So, typically, the neighbourhood of a feasible solution is specified by a small class of feasibility-preserving modifications/perturbations to its combinatorial structure. 2. A neighbourhood should be so defined as to ensure that it is polynomially searchable. Hence, unless we have a sophisticated non-exhaustive neighbourhood search algorithm, neighbourhoods should be small.
However, as we have argued in Sect. 1, this approach is not without issues. In particular, feasibility-preserving perturbations may not be easy to find, especially for NP-Hard-Feasible problems, while having small neighbourhoods implies a low likelihood of high-quality solutions.

Advanced Versions of Local Search.
A disadvantage of the basic version of local search is that it may stop at a local minimum of a poor quality, if too small a region of the feasible space is explored. To circumvent this outcome, advanced variants enabling an exploration of larger portions of the feasible space have been devised [7,29]. Those described here provide some mechanism to escape from the local minimum to "nearby" solutions and resume the search from there. They have been designed to accommodate situations, where local minima are not distributed uniformly in the feasibility space, but are rather clustered in close proximity [25].
The variable neighbourhood search approach uses multiple neighbourhoods to escape from local minima. It relies on the fact that a local minimum with respect to one neighbourhood need not be a local minimum with respect to another (if the latter is not contained in the former). The algorithm maintains a set of neighbourhood functions. Once a local minimum with respect to the current neighbourhood is reached, the neighbourhood is switched, and the search is resumed.
The hill climbing method allows the selection of a non-improving solution, once a local minimum is reached. Since the objective function no longer monotonically decreases, there is now a possibility of a cycle: a solution may be visited more than once forcing the search into an infinite loop. One can deal with this problem in various ways: ignore it and let the algorithm run until the timeout expires, use randomization, or introduce data structures that keep track of the search history and prohibit solutions that have already been encountered. The latter approach is referred to as tabu search.
Another idea is to use large neighbourhoods. This approach increases the size of the explored region and the likelihood of better solutions. However, large neighbourhood search may become intractable.
The iterated local search approach can be viewed as "a local search within a local search". In each iteration of the search, it uses a subsidiary search algorithm to explore iteratively a feasible sub-space. Once a local minimum is returned, a new search is initiated in a region, whose elements are obtained by "perturbing" the recent solution.
All the above approaches can be implemented within our LSSO framework. The key difference between LSSO and previous approaches is using SAT or Polosat as an oracle for both finding the initial solution and carrying out the neighbourhood search.

Local Search with SAT Oracle (LSSO)
This section introduces the main contribution of our paper. We propose using SAT as an oracle in local search algorithms to address the scalability and quality issues that arise in the classical local search algorithms, especially, given an NP-Hard-Feasible problem.
Given a combinatorial optimization problem, the first stage in designing an LSSO solution is expressing the problem as a COP.
In the second stage, the COP decision variables are translated to bit-vectors, and the feasibility constraints are translated to a BV formula (including any additional industrial requirements). One might experiment with several alternative formulations and select the one deemed best.
The third step is defining the so-called neighbourhood generators. A neighbourhood generator N (µ) accepts as an input a solution µ (that is, a model to the bit-vector formula, representing the COP), and generates neighbourhood constraints. The set of all the assignments which satisfy the feasibility and neighbourhood constraints constitutes the neighbourhood of the solution. Thus, finding such an assignment amounts to finding an element of the neighbourhood of µ.
A key ingredient of our methodology is the adoption of a neighbourhood concept, which differs significantly from the classical one, described in Sect. 2.4.2: 1. The neighbourhood need not be small and need not contain (only) elements "close" to the current solution.
2. Normally, N (µ) should generate constraints which ensure a cost lower than that of µ. If such a formulation is possible, then an iteration of the local search algorithm merely needs to find a model to these constraints in order to progress. 3. If the objective function is too complex to model in its entirety, a neighbourbood generator might attempt to ensure a better value for the objective function by imposing constraints on the objective function's sub-components. For example, when the objective function is a very large sum of bit-vector terms, one might impose constraints on the sum's terms or small partial sums thereof. 4. Notwithstanding the above, neighbourhood generators may support hill climbing, in which case, the constraints are so formulated as to admit non-improving solutions.
Note that, in our approach, neighbourhoods direct the search to "higher-quality" regions with respect to the current solution, regardless of the algorithmic difficulties of searching such regions. This is another key aspect of our approach: we trust SAT solvers to search complex sub-spaces efficiently.
Having discussed neighbourhoods, we are now ready to describe the simplest LSSO implementation: 1. A BV solver instance is created and the COP is provided to the solver. Specifically, we represent the COP's decision variables as bit-vectors, where the widths are chosen to accommodate the largest values. We provide the feasibility constraints to the solver as BV constraints. Then, we implement neighbourhood generators, which, given a feasible solution, return a set of BV constraints defining its neighbourhood. 2. The local search is carried out as follows: (a) The algorithm obtains an initial solution by asserting the feasibility constraints and asking the solver for a model. This model is set as the current solution µ. (b) The algorithm enters a loop, in which the solver operates in incremental mode.
In each iteration, the algorithm calls the neighbourhood generator with the current solution as input, to generate a list of BV constraints. These are provided to the solver, which is asked for a model. If a model α is found, µ is set to α. Otherwise, the algorithm terminates returning µ.
The neighbourhood constraints can be given to the solver as either assumptions or assertions. This leads to two types of search, providing a tradeoff between execution time and quality: 1. Non-speculative search: the neighbourhood constraints are passed to the solver as assertions. Once assertions are passed to the solver, they are enforced in all ensuing iterations. The search proceeds through a monotone sequence of decreasing neighbourhoods until a local minimum is reached. Thus, the search is localized and is relatively fast at the possible expense of quality. 2. Speculative search: the neighbourhood constraints are passed to the solver as assumptions. The neighbourhood constraints are valid only for one iteration. Thus, the current neighbourhood is not intersected with previous neighbourhoods and a larger portion of the feasibility space will be explored. The search is expected to be slower, since the SAT solver handles assumptions less efficiently than assertions [18], but the quality of resulting solution is expected to be better, since the search can explore a greater part of the feasibility space, especially so by variable neighbourhood search and hill climbing.
Alg. 1 depicts our implementation of LSSO. The algorithm receives four inputs. The Boolean inputs VN S, HC, and SPEC specify whether variable neighbourhood search, hill climbing, and speculative search are to be used. All combinations are possible, except that hill climbing requires speculative search. The input N max applies to variable neighbourhood search. It specifies an upper bound on the number of consecutive neighbourhood switches without finding a solution. If that bound is exceeded, the algorithm terminates with the current solution. To effect variable neighbourhood search, the algorithm uses a predefined list of neighbourhood generators N = [N 0 (µ), N 1 (µ) . . . ]. The first generator N 0 (µ) is considered the default and is used most of the time. The others are used to escape local minima.
Alg. 1 carries out iterated local search with Polosat as an oracle, where the observables are recommended to be set to the bits of the inputs of the objective function. One can also replace the Polosat invocation by an ordinary SAT invocation.

LSSO Algorithms for the Cell Placement Problem
This section presents our LSSO-based placement algorithms. All the algorithms are instantiations of Alg. 1 with different sets of parameters. The BV constraints are generated by translating the COP constraints, as discussed in Sect. 2.3. Each algorithm uses some of the neighbourhood generators defined in Sect. 4.1.
The algorithms are presented in Sect. 4.2. None of the algorithms define the target bit-vector explicitly, since they rely on local search instead of OBV solving. By default, the algorithms use Polosat as the oracle, where the observables comprise all the bits of the bit-vectors, representing the sizes of the nets, where the size of net n is given by the following bit-vector term (for every intermediate term and the resulting term n , its width is set to the minimal possible width which prevents an overflow, where the operators are zero-extended, whenever required): Let µ be a placement, that is, a model to the bit-vector formula representing the feasibility constraints. The neighbourhood N 1 (µ) is designed for a highly localized fast search at the possible expense of quality. To this end, the constraints corresponding to N 1 (µ) force a decrease of the objective function in a very constrained manner, so as to help the solver to come back quickly. N 1 (µ) consists of all of legal placements, for which all the nets are no bigger and at least one net is smaller than under µ, thus ensuring a lower cost. The constraints are: each net is no bigger n∈I ( n ≤ µ( n )) at least one net is smaller n∈I ( n < µ( n ) hill climbing constraints generator From now on, confine the search to the feasible space 2: solver.Assert(L) 3: current ← solver.Sat() find the initial solution 4: if ¬current then return None the problem is unsatisfiable Loop initialization 5: best ← current 6: stop ← ⊥ stopping condition 7: jump ← ⊥ indicates whether hill climbing should be attempted 8: i ← 0 current neighbourhood index 9: while ¬stop do Compute neighbourhood constraints 10: if HC ∧ jump then hill climbing is required 11: neighbourhood constraints := J (current) 12: else hill climbing is not required 13: If the mode is speculative, constraints are assumptions; otherwise they are assertions 14: if SPEC then 15: assertions solver.Assert(assertions) 19: next ← solver.P olosat(assumptions) 20: if next then found a solution 21: current ← next; i ← 0; jump ← ⊥ 22: if current.cost < best.cost then best ← current 23: continue Solution not found If we are in variable neighbourhood mode and the number of consecutive neighbourhood switches without a model has not exceeded the bound, move to next neighbourhood 24: if VN S ∧ (i < (Nmax − 1)) then 25: i ← i + 1 26: continue If we are in hill climbing mode, and have exhausted the bound on neighbourhood switches without getting a model, and hill climbing has not already been attempted in this iteration, attempt it in the next iteration 27: if HC ∧ ¬jump then 28: jump ← 29: continue If we got here, we are stuck and need to terminate 30: stop ← 31: return best 4.1.2 N 2 : a Family of Neighbourhood Generators. The N 2 family is designed for variable neighbourhood search. Each of its neighbourhoods strictly contains N 1 and allows the objective function to decrease in more ways. This implies higher quality solutions at the expense of slower convergence. To define the N 2 family, let α = I be the number of the nets and assume α ≥ 3. For each permutation σ of [1 . . . α] and positive number 2 ≤ d < α we define a neighbourhood function N 2 [σ, d](µ) as follows: Let n σ(1) , . . . n σ(α) be the permuted sequence of the nets. Partition this sequence into α/d segments of size d (last segment could be shorter). The neighbourhood N 2 [σ, d](µ) consists of all of legal placements, for which the sum of the net sizes of each segment is no bigger than under µ, and the sum of at least one segment is smaller. Note that this ensures a cost lower than the placement under µ. By choosing different pairs σ, d , one may obtain different neighbourhoods. The constraints are: Hill-climbing Neighbourhood Generator N 3 . N 3 is designed to implement hill climbing. We reason as follows: If the current placement is not a global minimum, there exists a placement with at least one smaller net. Hence, to tunnel away from the local minimum, we generate the following neighbourhood constraints: at least one net is smaller n∈I n < µ( n )

LSSO-based Algorithms for Placement
All the algorithms below are instantiations of Alg. 1; they use lists of neighbourhood generators, composed of the ones defined in Sect. 4.1, where hill climbing is carried out by using the neighbourhood generator N 3 . Due to project deadline constraints, we did not explore other combinations.

Experimental Results
We study the performance of the following algorithms within our placement tool: we chose many env spec hill clmb, since, as we shall soon see, it outperforms the other LSSO algorithms in a pairwise comparison). 3. virtual-best: represents the best result of the above algorithms per timeout.
We used an extensive set of 1200 proprietary industrial designs of various sizes and complexities. The sizes of the grids (where a grid size is the width N multiplied by the height M ) can be characterized as follows: a) Minimum size = 70; b) Maximum = 364000; c) Average ≈ 4643; d) Standard deviation ≈ 18829. We used machines with 32Gb of memory running Intel Xeon processors with 3Ghz CPU frequency.
We ran the algorithms for 600 seconds and measured the quality of the placement at different time intervals. Fig. 2 shows our main results. For each algorithm and time interval, Fig. 2 displays a score which represents the quality. The score is a real number between 0 and 1 inclusive, where the closer the score is to 1 the better. For each algorithm and time interval, the score is computed as follows: we compute the average value of the following score-per-instance: (the result of virtual-best in 600 sec.) / (the result of the current algorithm within the current time interval). Our conclusions: First, when using SAT as the oracle, Linear Search (ls no polosat) outperforms OBV-BS (bs no polosat), demonstrating that OBV-BS is not useful when the optimization target is a complex arithmetic expression (rather than a vector of lexicographically ordered bits, where each bit is a result of a separate calculation as in [19]). Based on this result, we preferred Linear Search over OBV-BS as the baseline algorithm.
Second, confirming the conclusion of [16], Polosat makes Linear Search substantially more efficient (compare ls to ls no polosat).
Third, and more importantly in the context of this work, our best novel LSSO algorithm even without Polosat (many env spec hill clmb no polosat) is almost as efficient as Linear Search with Polosat (ls), the latter being the state-ofthe-art in solving placement [16]. Moreover, the best Polosat-based LSSO algorithm (many env spec hill clmb) is significantly more efficient than both aforementioned algorithms. This result justifies the usage of both major components of our solution: LSSO-the high-level local search on top a satisfiability oracle, presented in this paper, and Polosat [16]-the low-level local search simulation with SAT.
Finally, the virtual best algorithm yields the absolutely best result, providing evidence that development of different LSSO algorithms pays off.
Additionally, Table 1 shows a pairwise comparison between our four Polosatbased LSSO algorithms. many env spec hill clmb outperforms the others. Table 2 offers a fine-grained comparison between our best novel LSSO algorithm many env spec hill clmb and the Polosat-based Local Search ls, the latter being the state-of-the-art in solving placement [16]. The comparison is provided per grid size category and for two different timeouts. LSSO improves the performance significantly for every input size category for both timeouts. Comparing the results for the two timeouts on the biggest instances shows that increasing the timeout makes the gap between LSSO and ls more significant, given large grids.
Finally, Table 3 shows the unique contribution of each algorithm to the virtual best in 600 sec. (we dismissed all the instances on which there was more than one bestperforming solver). Notably, each of the LSSO algorithms is a contributor. Surprisingly, many nbr nonspec contributes more than many env spec hill clmb, despite the latter algorithm outperforming the former in a pairwise comparison. A possible explanation is that we ran many nbr nonspec with Polosat only, while many env spec hill clmb was run twice with Polosat and SAT. Another surprising result is the significant contribution of many env spec hill clmb no polosat, second only to many nbr nonspec, implying that a SAT-based LSSO algorithm should be part of any parallel portfolio.  Table 1: Pairwise comparison between LSSO algorithms for the timeout of 600 sec. Each nonempty cell (r, c) contains a comparison between Algorithm R in row r and Algorithm C in column c. The value (w d l) in each non-empty cell is interpreted as follows: R outscored C on w instances; there was a draw on d instances; C outscored R on l instances.  Table  2: Comparing the best Polosat-based LSSO algorithm (many env spec hill clmb) to the Polosat-based Linear Search (ls), the latter comprising the previous state-of-the-art.

Conclusion
We have presented a new methodology for solving NP-hard combinatorial optimization problems, called Local Search with SAT Oracle (LSSO). Our approach can handle problems for which finding even one feasible solution is already NP-hard. LSSO applies local search which uses a SAT solver or the SAT-based optimization algorithm Polosat as an oracle. We have introduced a generic algorithm which integrates different local search schemes within the LSSO framework. Furthermore, we have implemented our approach in an industrial tool for solving the cell placement problem in VLSI and have shown that our new LSSO approach makes the tool substantially more efficient. Our tool has been successfully productized at Intel.

Algorithm
Contribution Algorithm Contribution many nbr nonspec 240 ls 33 many env spec hill clmb no polosat 181 many env spec 21 many env spec hill clmb 79 ls no polosat 12 single nbr nonspec 54 bs no polosat 8 Table 3: Unique contribution to the virtual best per algorithm (sorted by the contribution).


Compression-Based Methods of Time Series Forecasting †

: Time series forecasting is an important research topic with many practical applications. As shown earlier, the problems of lossless data compression and prediction are very similar mathemati-cally. In this article, we propose several forecasting methods based on real-world data compressors. We consider predicting univariate and multivariate data, describe how multiple data compressors can be combined into one forecasting method with automatic selection of the best algorithm for the input data. The developed forecasting techniques are not inferior to the known ones. We also propose a way to reduce the computation time of the combined method by using the so-called time-universal codes. To test the proposed techniques, we make predictions for real-world data such as sunspot numbers and some social indicators of Novosibirsk region, Russia. The results of our computations show that the described methods ﬁnd non-trivial regularities in data, and time universal codes can reduce the computation time without losing accuracy.

Introduction
The problem of time series forecasting is to estimate the future values of a process from a sequence of its observations. This task is important because it has many practical applications. Examples include predicting future stock prices and air temperature forecasting. Nowadays, there are many different approaches to solving this problem. Classical statistical models such as exponential smoothing and the autoregressive integrated moving average (ARIMA) model are very popular, highly accurate, and relatively easy to use. A detailed description of these methods can be found in [1,2]. Neural networks [3][4][5] are also widely used, especially on large datasets. However, there is no best method for all situations, and the development of new forecasting techniques remains relevant.
This work is based on an information-theoretic approach to time series forecasting. As is it was shown in [6], the problems of data compression and prediction are closely related and an asymptotically optimal method for predicting stationary stochastic processes can be based on a universal code (see also [7]). In [8] it was shown how any lossless data compression algorithm can be used to predict finite-alphabet and real-valued time series.
In this paper, we do not focus on asymptotic properties of algorithms and consider using this approach in practical situations. The main contributions of the work are summarized as follows. First, we describe how to use arbitrary data compressors for time series forecasting. Such compressors are a promising tool for forecasting because many of them are implementations of universal codes with numerous modifications to increase the level of compression in real-world situations. It is important to note that modern data compression techniques are based on several different approaches to universal coding: the Burrows-Wheeler Transform (BWT) [9], the Prediction by Partial Matching (PPM) [10] algorithm, the Lempel-Ziv family of algorithms [11,12], among others. Some algorithms [13][14][15] search for a compact context-free grammar that unambiguously represents the sequence for compression, one can see this approach as a kind of artificial intelligence. In previous works, only some theoretical possibilities of using universal codes for prediction without any practical applications were described [8], or the use of one universal code for prediction was presented [7].
Secondly, we propose an adaptive approach to time series forecasting, which is useful in situations when we do not know in advance which data compressor is optimal for a given time series.
Thirdly, we describe how the proposed techniques can be using to predict multivariate data. To test the proposed techniques, we make forecasts for real-world data such as sunspot numbers and some social indicators of Novosibirsk region, Russia.
It should be noted that our approach can be complemented with well-known time series transformation and adjustment techniques. For example, in this work, we use differencing and smoothing in the computational experiments.
The rest of the paper is structured as follows. In the next section, we describe the mathematical foundations of using arbitrary data compression techniques for forecasting finite-alphabet time series. After that, we describe generalizations of the model to realvalued and multivariate cases. Then we give some examples of the practical use of the presented techniques. Next, we propose an adaptive algorithm that can significantly reduce computation time when multiple data compressors are used. Further, we use the proposed algorithm to predict sunspot numbers. At the end of the paper, we describe the limitations of the proposed approach and make some conclusions.

Predicting Finite-Alphabet Time Series
Time series with finite alphabets are most convenient for forecasting using data compression algorithms. Suppose we have a sequence X = x 1 , x 2 , . . . , x t , x i ∈ A, where A is a finite set (an alphabet), and we want to give a prediction for x t+1 , x t+2 , . . . , x t+h , h ∈ Z + .
Denote as A n the set of all sequences of lengths n over A and A * = ∞ i=0 A i . A uniquely decodable data compression method (or code) ϕ is a set of mappings ϕ n : A n → {0, 1} * , n = 1, 2, . . ., such that for any sequence of words x 1 , x 2 , . . . , x m , x i ∈ A n , m ≥ 1, the sequence ϕ n (x 1 ), ϕ n (x 2 ), . . . , ϕ n (x m ) can be uniquely decoded as x 1 , x 2 , . . . , x m . The compressed size (in bits) of sequence α we denote as |ϕ(α)|. We can get a probability distribution on A n using ϕ by (1) A code ϕ is called universal if for any stationary and ergodic measure P with probability 1, and where H(P) is the entropy rate [16] of P, E( f ) is the expected value of f . In [8] it was shown that if P is a stationary and ergodic stochastic process and ϕ is a universal code, (1) in certain sense is a nonparametric estimate of the unknown probability measure P(X). More precisely, the following theorem was proved.
Theorem 1. If P is a stationary ergodic measure and ϕ is a universal code, then the following equalities hold: . . , x t ))) = 0 with probability 1,
Suppose that A is a finite set of integers, P ϕ (x t+j = a|X) = ∑ Y=(y 1 ,...,y j−1 ,a,y j+1 ,...,y h )∈A h P ϕ (Y|X). We can give a point forecast asx t+j = ∑ a∈A aP ϕ (x t+j = a|X) (i.e., compute the mean over the marginal distribution for step j).
Example 1. Consider the following sequence: Suppose we want to make a forecast for its next two values using gzip [17] as ϕ. Assume that the alphabet A of the underlying process is {0, 1}. We need to compress every sequence of the form XZ, where Z ∈ A 2 . The compression results along with the calculated conditional probabilities P gzip (Z|X) are presented in Table 1.
Suppose we have a finite set of data compression algorithms F = {ϕ 1 , ϕ 2 , . . . , ϕ k }. If we do not know which ϕ ∈ F is the best predictor for a given series X, we can mix the conditional probability distributions yielded by each compressor from F by where Note that (3) works in such a way that the highest probability gets Y ∈ A h such that |ϕ s (XY)| = min Z∈A h ,ϕ∈F |ϕ(XZ)| for some ϕ s ∈ F, i.e., (3) selects the best compressor ϕ s "automatically".

Predicting Real-Valued Time Series
Many time series that can be found in practice are sequences of real numbers. To predict such a series using data compression algorithms, we need to convert it to a sequence of integers (with loss of information, obviously). This process of conversion is known as quantization. Consider a sequence X = x 1 , x 2 , . . . , x t , where x i ∈ R. Denote its minimal and maximal elements as m and M respectively: Probably the simplest way of conversion is to split [m; M] into a finite number n of disjoint numbered intervals {q 1 , q 2 , . . . , q n } of equal length and replace each x i with its corresponding interval number: if x i ∈ q j , then replace x i with j. Later, we can perform the inverse conversion, replacing the indices, for example, with the medians of the corresponding intervals. Now, we consider the question how to select the number of intervals n. On the one hand, if this value is too small, some important regularities may be missing in the converted series. On the other hand, if n is too large, a data compressor may not be able to capture the regularities due to noise in the data. One possible solution of this problem is to employ an approach similar to that used in (3). Let n be some positive power of 2: l = log 2 n is an integer greater than zero. We can mix the probability distributions, obtained using partitions into 2 k intervals, 1 ≤ k ≤ l, with some weights. The partition yielded the smallest code length will have the greatest impact on the final result. Denote the number of interval that contains x i in partition into 2 k intervals as x [k] i . Then the mixed probability distribution can be defined as where A k = {0, 1, . . . , 2 k − 1} is a set of interval numbers, ω k -non-negative weights, l ∑ k=1 ω k = 1.

Predicting Multivariate Data
In some practical situations, it is required to predict several related series. In such cases, we can try to utilize the connection between them to improve the accuracy of our forecasts. For example, air temperature and barometric pressure are related, and it might make sense to predict them together. Many univariate time series forecasting techniques have generalizations to the multivariate (or vector) case [18][19][20][21].
The approach based on data compression algorithms can be used in a multivariate setting too. Suppose we have a vector time seriesX =x 1 , Let's find the minimal and maximal values for each coordinate: As in the previous section, we split each interval [m j ; M j ] into a finite number of intervals n and thus obtain n d d-cubes. Then we need to number these cubes and replace each pointx i with the number of the cube it falls into. As a result, we get an integer sequence that can be predicted using the previously described methods.

Experiments
In this section, we use the proposed method to predict real-world data. Along with the point forecasts (which are the means of the corresponding distributions), we provide 95% confidence intervals for them. To estimate these intervals, we used the following strategy. Suppose that we want to make an h-step ahead prediction for time series X = x 1 , x 2 , . . . , x t . Let s be t/2 . For each series X s , X s+1 , . . . , X t−h , where X i = x 1 , x 2 , . . . , x i , we made an h-step forecastX h i =x i+1 ,x i+2 , . . . ,x i+h and calculated the residuals for each step: Then we computed the standard deviations of the residuals and we want to make a 2-step forecast. As was explained previously, we make 2-step ahead predictions for the series X 2 = 1.1, 1.2 and X 3 = 1.1, 1.2, 1.3.
We mixed these compressors with the same weights, that is, in (3) ω i = 1/8. Several simple techniques of data transformation also were used. First, to remove trends in data we took the first difference: y i = x i − x i−1 . Secondly, we used smoothing: Let us move on to forecasting some indicators of Novosibirsk region. In the first example, we predict the annual number of unemployed persons in the region. In Figure 1a this series along with our 4-step forecast with confidence intervals is presented. Since this series is a real-valued one, we used partitions of its range of values into 2, 4 and 8 intervals in (4). According to our forecast, in the next four years the indicator will be higher than the current level. In the second example, we predict the average annual monetary income of the population in the region. The results are shown in Figure 1b. The forecast shows that in the next 4 years the trend will continue and by 2023 the indicator will be near 35,000 rubles per year.

Adaptive Method of Forecasting
As noted in Section 2.1, multiple data compression algorithms can be combined into one method of forecasting (almost) without loss of accuracy using (3). The only issue is computation time-we need to compress all sequences with every data compression algorithm that we include in our combination. In this section, we consider a way to significantly reduce computation time while maintaining near-optimal accuracy. This can be achieved using the so-called time-universal codes.

Time-Universal Codes
Suppose we want to compress a sequence X = x 1 , x 2 , . . . , x n , x i belongs to some finite alphabet A. Also suppose we have a finite set of data compression algorithms F = {ϕ 1 , ϕ 2 , . . . , ϕ k } and we want to compress X with ϕ s ∈ F that yields the smallest code length. The obvious way to find ϕ s is to compress X with each ϕ i and then to choose the best one. After that, we need to store s along with ϕ s (X) in order to be able to decompress the original data. Binary representation < q > of any integer q from 0 to k − 1 requires no more than log 2 k bits, therefore our final code word would be < s > ϕ s (X) with length log 2 k + |ϕ s (X)| bits. This approach works well but requires additional time proportional to k. The goal is to compress X with the best compressor using relatively small additional time.
Time-universal codes were proposed in [25] and allow making the portion of extra time asymptotically as small as desired. Let v i be the time ϕ i spends on encoding a single letter of X, v = max 1≤i≤k v i . Thus, an upper bound for time needed to encode X using any ϕ ∈ F is T = vn. Denote as δT the amount of extra time we have to select a close to optimal data compressor from F for X, where δ is a positive constant. The total time of selection and compression becomes no more than T + δT = T(1 + δ). Time-adaptive and time-universal codes are defined as follows.
Definition 1. Any method of data compression that encodes a sequence x 1 , x 2 , . . . , x n , n > 0, x i ∈ A by a binary word of the length log 2 k + |ϕ s (x 1 , x 2 , . . . , x n )| for some ϕ s ∈ F and the time of encoding is not greater than T(1 + δ) is called as time-adaptive code and denoted asΦ δ compr .

Definition 2. If for a time-adaptive codeΦ δ compr the following equation is valid
this code is called time-universal.
In [25] was proposed the following simple algorithm and proved that it yields a time universal code:
In this work, we implemented this algorithm and used it to predict real-world time series.

Experiments
In this section, we consider sunspot numbers forecasting. Sunspots are temporary spots on the Sun with reduced surface temperature that appear darker than the surrounding areas. The sunspot number time series is provided by the WDC-SILSO (World Data Center-Sunspot Index and Long-term Solar Observations) [26]. The daily, monthly, yearly and 13-month smoothed monthly total sunspot numbers can be found on the SILSO site http://www.sidc.be/silso/datafiles. Here we used the monthly mean total sunspot number series.
We accessed the SILSO site on 20 June 2020, at that time the series had 3257 observations. The entire series is shown in Figure 2. It has an obvious cyclical component and it's known that the approximate cycle lengths is 11 years [27].
The Space Weather Services (SWS) of the Australian Bureau of Meteorology issues forecasts for the WDC-SILSO sunspot numbers and maintains an archive of this forecasts (http://listserver.ips.gov.au/pipermail/ips-ssn-predictions/), so we can compare our forecasts with them. In July 2015, SILSO adjusted the original data and SWS moved to the new version in the forecasts in February 2016. We started making predictions from that date, thus after each element with the number greater than 3205 (which corresponds to January 2016), before adding it to the end of the series, we made a 4-step forecast. In order to select a (close to) optimal compressor for this series, we tried to use from 10% to 100% of its values with step 10%. For instance, when 10% of the values were used, δ was 8× 0.1 = 0.8 since we had 8 algorithms. In Section 2.4, the maximal number of intervals we considered when quantizing was 8; in this section, we increased this value to 16. In everything else, we used the same methodology as in Section 2.4.
For brevity, let us denote the combination of all 8 algorithms, obtained using (3), as joint method.Code lengths for different compressors, obtained by compressing the entire series (100% of values), are presented in Table 2. We can see that by code length zstd is the best compressor for this series and hence it will have the greatest contribution to the result of the joint method. For 10% trough 40% of the series values the best compressor by code lengths is ppmd, then for 50-100% zstd becomes the best compressor. The code length for ppmd and zstd are shown in Table 3. To assess the accuracy of our forecasts, we calculated the mean absolute error (MAE), defined as the mean of the absolute values of the residuals, for each step 1-4. The MAEs of ppmd's and zstd's forecasts are presented in Table 4. We can see that ppmd was a bit more accurate than zstd. We tried to add more values for prediction (from 2800 to 3205), and in that case zstd was slightly more accurate for steps 1-3 (the mean absolute errors are presented in Table 5). For step 4, ppmd was more accurate. We ran the program 5 times using all 8 algorithms and 5 more times using the adaptive algorithm (50% of the values were used to select the best compressor). The average time required to compute one 4-step forecast was 9945.46 s. in the former case and 564.63 s. in the latter. Thus, we reduced the computation time by more than 17 times without loss of accuracy.

Model Limitations
The presented model has some limitations. First, it is suitable for prediction of data with no trend pattern. But, as was shown in Section 2.4, to overcome this limitation some additional techniques such as differencing or time series decomposition [28] can be used. Secondly, the computational complexity of the model is high: if we make forecast for h steps ahead, |A| h sequences have to be compressed. This means that it cannot be directly applied in long-term forecasting or in a setting when it is required to compute predictions very quickly. While sacrificing accuracy, we can still make long-term forecasts using the following approach. Suppose we have a series X = x 1 , x 2 , . . . , x t and we want to predict its next h values. For simplicity, assume that t and h are even numbers, but the generalization is obvious. We can split X into two separate time series X odd = x 1 , x 3 , . . . , x t−1 and X even = x 2 , x 4 , . . . , x t . Then we predict x t+1 , x t+3 , . . . , x t+h−1 using X odd and x t+2 , x t+4 , . . . , x t+h using X even . This allows us reduce the number of sequences to compress from |A| h to 2|A| h/2 . It is clear that we can go further and split X into more than two series. Another obvious way to speed up computations is to choose small n when quantizing.

Discussion
In our opinion, the results of computations show that the proposed method has good accuracy. The adaptive approach to forecasting can significantly reduce computation time without loss in effectiveness if multiple compressors are used.
It is important to note that some time series can contain complex regularities. For example, the dynamics of financial time series can be influenced by the participants of the corresponding processes, and this can lead to the emergence of subtle patterns in the data. Another example is the time series arising in the study of space objects. Such series can contain various nonstationarities like attenuation and non-periodic changes. Thus, in our opinion, forecasting methods that allow finding unusual patterns may be of interest to many researchers. Patterns of the indicated types that occur in text data can be found by modern archivers. Therefore, we believe that the use of data compressors implicitly allows the use of various techniques of finding patterns, including artificial intelligence methods beyond neural networks.
In further research, a more elaborate strategy for optimal algorithms selection in the adaptive method can be developed. For example, one can use multidimensional optimization techniques for this purpose.

Conclusions
It turns out that from a mathematical point of view, data compression and prediction can be thought of together, allowing ideas from one area to be used in another. From a practical point of view, many of the data compression techniques implemented in software can be used in time series forecasting. Efficient data compression algorithms to be developed in the future can be combined with existing ones using the described approach.  Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://novosibstat.gks.ru/folder/31734 (average number of unemployed persons in Novosibirsk region), https://novosibstat.gks.ru/folder/31847 (average annual monetary income of the population in Novosibirsk region), http://www.sidc.be/silso/datafiles (monthly mean total sunspot numbers), http://listserver.ips.gov.au/pipermail/ips-ssn-predictions/ (SWS observed and predicted solar indices).

Conflicts of Interest:
The authors declare no conflict of interest.