EDA T9/T10 — Exhaustive Search and Generation

Ilario Bonacina
ilario.bonacina+EDA@upc.edu

This presentation is an adaptation of material from Antoni Lozano and Conrado Martinez

Agenda

Last theory class

  • Algorithms on Graphs
    • BFS
    • Dijkstra
    • Minimum Spanning Trees

Today

  • Exhaustive Search and Generation
    • Pruning
    • Branch & Bound
    • Several examples

Before we start: Questions on previous material?

Intro

In many computational problems we must exhaustively explore a solution space.

Many times it is the only way we know to solve the problem at hand.

For instance,

  • Generating all combinatorial objects of a given size
    • bitstrings
    • permutations
    • trees
    • graphs
    • …
  • Finding a non-menacing configuration of n queens in a n\times n chessboard (n-Queens)
  • Finding a shortest Hamiltonian cycle in a weighted graph (Traveling Salesman Problem)
  • …and many more

Generation

Example: Strings of Zeros and Ones

Input:
a natural number n

Output:
all sequences of 0s and 1s of length exactly n



What is the size of the output?

#include<iostream>
#include<vector>

using namespace std;

void binary(vector<int>& A, int k) {
    if (k == A.size()){ // base case
      for(int j = 0; j < k; ++j) cout << A[j];
      cout << endl; 
      return;
    }  
    // inductive case
    A[k] = 0; 
    binary(A,k+1);
    A[k] = 1; 
    binary(A,k+1);
  }
  
  int main(){
    int n;
    cin >> n;
    vector<int> A(n);
    binary(A,0);
  }
1
i is the next position in A we will assign

A first blueprint for backtracking

void backtracking(int k, ...other parameters...){
    if (k == base_case){// base case of the recursion
        ...
        return;
    }
    // inductive case
      ...
}
1
Here we arrived to a full solution. Do here whatever you are supposed to do with a full solution, e.g. write it to cout.
2
Assume you already have a partial solution, for instance a vector whose entries up to k-1 represent a partial solution.

We want to extend the partial solution to a partial solution up to k, in all the possible allowed ways to extend the partial solution to a partial solution one element larger. (Sometimes this exploration is only on 2 elements, e.g. true or false and it is not done explicitly with a for loop)

Example: Zeros and Ones with a fixed number of ones

Input:
two natural numbers n and \ell with \ell\leq n

Output:
all sequences of 0s and 1s of length n with exactly \ell ones

What is the size of the ouput?

Answer \binom{n}{\ell} which is upper and lower bounded as follows \left(\frac{n}{\ell}\right)^\ell\leq \binom{n}{\ell}\leq \left(\frac{en}{\ell}\right)^\ell Hence for instance when \ell=\frac{n}{2}, \binom{n}{n/2} is exponentially large, unlike when say \ell=1.

A first solution

#include<iostream>
#include<vector>

using namespace std;

void strings(vector<int>& A, int k, int l) {
  if (k == A.size()) {
    int c = 0;
    for (int x : A) c += x; // Counts 1s
    if (c == l){
      for(int j = 0; j < k; ++j) cout << A[j];
      cout << endl;
    }  
    return;
  }
  A[k] = 0; strings(A, k+1, l);
  A[k] = 1; strings(A, k+1, l); 
}

int main(){
  int n, l; 
  cin >> n >> l;
  vector<int> A(n);
  strings(A,0,l);
}
1
A: partial string (size n)
k: first non-filled cell in A
l: total number of 1s we want

Do we really need to count from scratch the number of 1’s in A every time?

  • We could keep, at every moment, the number of 1s of the partial string.

  • This forces us to use one more parameter in the procedure.

#include<iostream>
#include<vector>

using namespace std;

void strings2(vector<int>& A, int k, int u, int l) {
  if (k == A.size()) {
    if (u == l){
      for(int j = 0; j < k; ++j) cout << A[j];
      cout << endl;
    }
    return;
  }
  A[k] = 0; strings2(A, k+1, u, l);
  A[k] = 1; strings2(A, k+1, u+1, l);
}
int main(){
  int n, l; 
  cin >> n >> l;
  vector<int> A(n);
  strings2(A,0,0,l); 
}
1
A: partial string (size n)
k: first non-filled cell in A
u: number of 1s in A[0...k-1] (already filled)
l: total number of 1s we want

Can we detect situations in which a partial string cannot be extended to a total string with exactly \ell 1s?

If we can do that we can avoid to visit branches in the search tree, i.e. we prune the search space.

#include<iostream>
#include<vector>

using namespace std;

void strings3(vector<int>& A, int k, int z, int u, int l) {
  if (k == A.size()){
    for(int j = 0; j < k; ++j) cout << A[j];
    cout << endl;
    return;
  }  
  if (z < A.size() - l) {
    A[k] = 0; strings3(A, k+1, z+1, u, l); 
  }
  if (u < l) {
    A[k] = 1; strings3(A, k+1, z, u+1, l); 
  }
}
int main(){
  int n, l; 
  cin >> n >> l;
  vector<int> A(n);
  strings3(A,0,0,0,l); 
}
1
A: partial string (size n)
k: first non-filled cell in A
u: number of 1s in A[0...k-1] (already placed)
z: number of 0s in A[0...k-1] (already placed)
l: total number of 1s we want
2
not all 0s placed
3
not all 1s placed

Example: Permutations

Input:
a natural number n

Output:
all possible permutations of n


What is the size of the output?
Answer n! 

#include<iostream>
#include<vector>

using namespace std;

void write_perm(int n, vector<int>& A, int k, vector<bool>& used) {
    if (k == A.size()){
        for(int j = 0; j < k; ++j) cout << A[j];
        cout << endl;
        return;
    }  
    for (int i = 1; i <= n; ++i) {
      if (not used[i]) {
        A[k] = i;
        used[i] = true;
        write_perm(n,A,k+1,used);
        used[i] = false;
      }
    }
  }
  int main() {
    int n; 
    cin >> n;
    vector<int> A(n);
    vector<bool> used(n+1,false);
    write_perm(n,A,0,used); 
  }
1
What happens if we don’t have this line of code?

A generic backtracking algorithm

  • The solution space is organized as a configuration tree. Each node or tree configuration is represented with a vector A = (a_1, a_2, . . . , a_k ) that contains the choices already made (the partial solution).

  • Vector A grows in the forward phase by choosing a_{k +1} from a set S_{k +1} of candidates (exploration in depth).

  • A is reduced in the backtracking phase (backtrack).

A second blueprint for backtracking

void backtracking(int k, ...other parameters...){
    if (k == base_case){// base case of the recursion
        ...
        return;
    }
    // inductive case
    ...for loop exploring all possible ways of extending...
}

Inside the for loop exploring all the possible extensions, check whether the way you are trying to extend is allowed (this really depends on the problem at hand).

If it is allowed:

  • Extend the partial solution.
  • Call backtracking with the the auxiliary structures upadted, i.e. call backtracking(k+1, ...other parameters updated...);
  • Undo the update on the auxiliary structures to allow trying other branches of execution

Common errors

  • Trying to do more than one check in the if condition of the base case. Or forgetting to have a return; at the end. This will likely result in some unwanted execution, for instance entering in the inductive case when k > base case and perhaps going out the range of a vector.

  • Assuming that the base case is k=0. This is not true!
    The base case of the recursion is when we have a full solution. k=0 is signaling the opposite, the partial solution is empty, we have not made yet any progress towards a full solution.

  • In the inductive case, trying to update k with ++k or updating other quantities in the same way instead of calling backtracking(k+1, ...).

Common errors

  • In the inductive case, not assuming you already have a partial solution up to k-1, and therefore trying to compute it again, typically with some for loop.

  • A common error is forgetting to undo the marking after the call to backtracking(k+1, ...). This will result in exploring only one possible branch (one way of extending partial solutions). This might lead or not to a full solution.

Example: n-Queens

Input:
a positive integer n

Output:
write all possible positions of n queens in an n\times n chessboard such that no queen attack another one (i.e. each row, column or diagonal contains at most one queen).

A First solution

  • finds all solutions using backtracking
  • extends the partial solution whenever it is legal (can be extended to a complete solution)

Worst-cast cost in time: \Theta(n^n)

Use vector<int> t; to store the positions of the queens:

The queen in row i is in column t[i].

When two queens can attack each other?

The i-th queen and the j-th queen can attack each other if

  • t[i] == t[j] (same column)
  • t[i]-t[j] == i-j (same diagonal)
  • t[i]-t[j] == j-i (same anti-diagonal)
bool attacked(int i) {
  for (int k = 0; k < i; ++k)
    if (t[k] == t[i] or t[k] - k == t[i] - i or t[k] + k == t[i] + i)
      return true;
  return false;
}
1
This function checks given the positions of the queens up to i, if the i-th queen is attacked by the previous ones.

To write the output we use write()

void write() {
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j)
      cout << (t[i] == j ? "Q" : ".");
      cout << endl;
    }
  cout << endl;
}
Q.......
....Q...
.......Q
.....Q..
..Q.....
......Q.
.Q......
...Q....

Using the functions write and attacked then we can output all the solutions as follow:

void queens(int i) {
  if (i == n) write();
  else
    for (int j = 0; j < n; ++j){
      t[i] = j;
      if (not attacked(i))
        queens(i+1);
    }
}
1
j is the column for queen of row i

A Second solution: using marking

Instead of computing whether a position is attacked every time we can use a marking strategy remembering which rows/columns/diagonals are not allowed since they already have a queen in it.

  • mark used columns (vector<bool> of size n)
  • 2n-1 anti-diagonals, identified by i+j. I.e. numbers between 0 and 2n-2. Use a vector<bool> of size 2n-2.
  • 2n-1 diagonals, identified by i-j. This gives numbers between -(n-1) and n-1. It is more convenient to use i-j+(n-1), which gives numbers between 0 and 2n-2. Use a vector<bool> of size 2n-2.

#include <iostream>
#include <vector>
using namespace std;
int n;
vector<int> t;
// mc[j] == queen at column j
// md1[k] == queen at diagonal i+j = k, etc.
vector<int> mc, md1, md2;

void queens(int i);

int main() {
  cin >> n;
  t = vector<int>(n);
  mc = vector<int>(n, false);
  md1 = md2 = vector<int>(2*n-1, false);
  queens(0);
}

queens function

void queens(int i) {
  if (i == n){
    write();
    return;
  }

  for (int j = 0; j < n; ++j) // j is column for queen of row i
    if (not mc[j] and not md1[i+j] and not md2[i-j + n-1]) {
      t[i] = j;
      mc[j] = md1[i+j] = md2[i-j + n-1] = true;
      queens(i+1);
      mc[j] = md1[i+j] = md2[i-j + n-1] = false;
  }
}

If we only wanted to print one solution?

bool queens(int i) {
  if (i == n){
    write();
    return true;
  }

  for (int j = 0; j < n; ++j) // j is column for queen of row i
    if (not mc[j] and not md1[i+j] and not md2[i-j + n-1]) {
      t[i] = j;
      mc[j] = md1[i+j] = md2[i-j + n-1] = true;
      if (queens(i+1)) return true;
      mc[j] = md1[i+j] = md2[i-j + n-1] = false;
  }
  return false;
}

If we wanted to count solutions?

int queens(int i) {
  if (i == n){

    return 1;
  }
  int res = 0;
  for (int j = 0; j < n; ++j) // j is column for queen of row i
    if (not mc[j] and not md1[i+j] and not md2[i-j + n-1]) {
      t[i] = j;
      mc[j] = md1[i+j] = md2[i-j + n-1] = true;
      res += queens(i+1);
      mc[j] = md1[i+j] = md2[i-j + n-1] = false;
  }
  return res;
}

Example: Latin Squares

Input:
a positive integer n

Output:
all Latin squares of size n, i.e. all n\times n grids filled with n distinct symbols and each symbol appear once in every row and column.

Backtracking solution with markings

#include <iostream>
#include <vector>
using namespace std;
int n;
vector<vector<int>> q;
vector<vector<bool>> r;
vector<vector<bool>> c;

void latin_square(int i, int j);

int main () {
  cin >> n;
  q = vector<vector<int>>(n, vector<int>(n));
  r = c = vector<vector<bool>>(n, vector<bool>(n, false));
  latin_square(0, 0);
}
1
q[i][j] equals value at row i, column j
2
r[i][v] equals whether row i already uses value v
3
c[j][v] equals whether column j already uses value v

Cost: \mathcal O(n^{n^2} )

latin_square function

// Find all Latin squares completing from (i, j)
void latin_square(int i, int j) {
  if (i == n) 
    return write();
  if (j == n) 
    return latin_square(i+1, 0);
  for (int v = 0; v < n; ++v) {
    if (not r[i][v] and not c[j][v]) {
      r[i][v] = c[j][v] = true;
      q[i][j] = v;
      latin_square(i, j+1);
      r[i][v] = c[j][v] = false;
    }
  }
}
void write() {
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      cout << q[i][j] << ’\t’;
    }
    cout << endl;
  }
  cout << endl;
}

More search problems

Example: Knight’s tour

Input:
a positive integer n and i,j\leq n

Output:
a sequence of positions of the n\times n board s.t.

  • they start from position (i,j)
  • they are covering all the positions of the board and
  • any two consecutive positions are reachable using a knight’s move.

A backtracking solution

#include <iostream>
#include <vector>
using namespace std;

int n;
vector<vector<int>> t;

bool possible(int i, int j, int s);

int main() {
  int i, j;
  cin >> n >> i >> j;
  t = vector<vector<int>>(n, vector<int>(n, -1));
  cout << possible(i, j, 0) << endl;
}
1
t[i][j] == k if in the k-th jump we arrive at (i,j)
-1 if we have not arrived at (i,j) yet
2
Can we fill board starting at (i, j) having done s jumps? 

vector<int> di = {1, 1, -1, -1, 2, 2, -2, -2};
vector<int> dj = {2, -2, 2, -2, 1, -1, 1, -1};

bool possible(int i, int j, int s) {
  if (i >= 0 and i < n and j >= 0 and j < n and t[i][j] == -1){
    t[i][j] = s;
    if (s == n*n-1) return true;
    for (int k = 0; k < 8; ++k)
      if (possible(i + di[k], j + dj[k], s+1)) return true;
    t[i][j] = -1;
  }
  return false;
}

Example: Hamiltonian Cycle

Input:
a graph G

Output:
True/False according whether G has a Hamiltonian cycle, i.e. a cycle that visits each vertex of G exactly once


vector<vector<int>> G;
int n; //number of vertices
vector<int> s;
bool found;
vector<int> sol;

void hamilton(int v, int t);

int main(){
  ...read the input graph G...
  n = G.size();
  s = vector<int>(n, -1);
  found = false;
  hamilton(0,1);
  cout << found << endl;
}
0
The input graph. Assume it is connected and represented via adjacency lists
1
s is the partial solution we build. s[i] is the vertex after i in the partial solution (−1 if not used yet)
2
whether we have already found a cycle
3
a hamiltonian cycle (if found)
4
v = last vertex in the path
t = number of vertices in the path 

void hamilton(int v, int t) {
  if (t == n) {
    if (G[v][0] == 0) {
      s[v] = 0;
      found = true;
      sol = s;
      s[v] = -1;
    }
    return;
  }
  for (int u : G[v]) {
    if (s[u] == -1) {
      s[v] = u;
      hamilton(u, t+1);
      s[v] = -1;
      if (found) return;
    } 
  } 
}
1
should make sure path s can be closed to a cycle

Optimization problems

Example: Knapsack

Input:
a positive integer C and n pairs of numbers (v_1,w_1),\dots,(v_n,w_n)

Output:
a subset S of \{1,\dots,n\} such that

  • \sum_{i\in S} w_i\leq C and
  • \sum_{i\in S}v_i is maximum, i.e. for every other subset T such that \sum_{i\in T} w_i\leq C we have \sum_{i\in S}v_i\geq \sum_{i\in T}v_i

First solution: prune when capacity is exceeded

#include <iostream>
#include <vector>
using namespace std;

int c; // Capacity
int n; // Number of objects
vector<int> w; // Weights
vector<int> v; // Values
vector<int> s; // (Partial) Solution we are building
// s[i] == 1 iff i-th object is chosen
int bv = -1; // Best value so far
vector<int> bs; // Best solution so far

void opt(int k, int swp, int svp);

int main() {
  cin >> c >> n;
  w = v = s = vector<int>(n);
  for (int& x : p) cin >> x;
  for (int& x : v) cin >> x;
  opt(0, 0, 0);
  cout << bv << endl; 
}

void opt(int k, int swp, int svp) { // swp: sum weights partial
  if (swp > c) return; // Exceed capacity: do not continue
  if (k == n) {
    if (svp > bv) { // Improve best solution so far
      bs = s;
      bv = svp;
    }
    return;
  }
   // Discard obj. k
  s[k] = 0; 
  opt(k+1, swp, svp );
  // Choose obj. k
  s[k] = 1; 
  opt(k+1, swp + w[k], svp + v[k]); 
}

Second solution: Branch and Bound

Prune also when we cannot improve the best cost found so far

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

...same global parameters as before...

// svr: sum values remaining
void opt(int k, int swp, int svp, int svr);

int main() {
cin >> c >> n;
p = v = s = vector<int>(n);
for (int& x : p) cin >> x;
for (int& x : v) cin >> x;
sum = accumulate(v.begin(), v.end(), 0);
opt(0, 0, 0, sum);
cout << bv << endl;
}

// svr: sum values remaining
void opt(int k, int swp, int svp, int svr) {
  if (swp > c or svp + svr <= bv) return;
  if (k == n) {
    bs = s;
    bv = svp;
    return;
  }
  s[k] = 0; 
  opt(k+1, swp, svp, svr - v[k]);
  
  s[k] = 1; 
  opt(k+1, swp + w[k], svp + v[k], svr - v[k]);
}

Example: Travelling Salesman Problem (TSP)

The Travelling Salesman Problem (TSP) consists in, given a network of cities, finding the order in which we should visit the cities such that:

  • we start and finish in the same city
  • we visit each city exactly once, and
  • the total travelled distance is as small as possible

  • It is one of the most famous combinatorial optimization problems
  • It is important from the theoretical point of view
  • It has practical applications in
    • planning
    • logistics
    • microchips manufacturing
    • DNA sequence
    • astronomy
    • …

class TSP {
int n; matrix<int> M; vector<int> s; vector<int> best_sol; double best_cost;

void recursive (int v, int t, double c) {
  if (t == n) {
    c += M[v][0];
    if (c < best_cost) {
      best_cost = c;
      best_sol = s;
      best_sol[v] = 0;
    }
  return;  
  }
  for(int u = 0; u < n; ++u)
    if (u != v and s[u] == -1) {
      if (c + M[v][u] < best_cost) {
        s[v] = u;
        recursive(u, t+1, c + M[v][u]);
        s[v] = -1;
      } 
    } 
}
public:
  TSP (matrix& M) {
    this->M = M;
    n = M.rows();
    s = vector<int>(n, -1);
    best_sol = vector<int>(n);
    best_cost = infinite;
    recursive(0, 1, 0);
  }
  vector<int> solution ( ) { return best_sol; }
  int next (int x) { return best_sol[x]; }
  double cost ( ) { return best_cost; }
};
1
n: number of cities
M: matrix of distances
s is the partial solution we build. s[i] is next city after city i (−1 if not yet used)
best_sol: best solution found so far
best_cost: cost of best solution found so far
2
v : last vertex in partial solution s
t : size of the path given by s
c : cost of s

Current best solution pruning (CBSP)

Current best solution pruning (CBSP) is a technique that can be applied whenever we want to solve optimization problems (maximization or minimization). Say, our goal is to find a solution of minimum cost.

Let x be the current solution, and best_cost the cost of the best solution found so far. We compute an estimation of the cost of the best solution that we can find from x (estimated_cost). It must be a lower bound the cost of the best solution that we could find from x (real_best_cost).

Then if estimated_cost > best_cost we can prune at x and avoid exploring its descendants in the recursion tree: no solution found there will be better than the one we already have.

Upcoming classes

Wednesday
Lab class P5 on exhaustive search (do the exercises before class).
Next Monday
Theory class T11 on Intractability




Questions?