# jemdoc: addcss{rbh.css}, addcss{jacob.css}
= hard problems

~~~
{}{raw}
 <a href="#bad">bad</a> &nbsp;
 <a href="#hard">hard</a> &nbsp;
 <a href="#knapbf">knap:bf</a> &nbsp;
 <a href="#knapdp">knap:dp</a> &nbsp;
 <a href="#knapdpv">knap:dpv</a> &nbsp;
 <a href="#setcover">set cover</a> &nbsp;
 <a href="#sat">sat</a> &nbsp;
~~~

#~~~
#{}{raw}
#<h2 id="bad">baby, you with the bad boys now</h2>
#~~~
#
#~~~
#from Out of Sight by Elmore Leonard
#~~~

~~~
{}{raw}
<h2 id="hard">hard problems</h2>
~~~

~~~
for many problems, no known polytime algorithm
- [https://en.wikipedia.org/wiki/Graph_isomorphism graph isomorphism]
- [https://en.wikipedia.org/wiki/Knapsack_problem 0-1 knapsack]
- [https://en.wikipedia.org/wiki/Set_cover_problem set cover]
- [https://en.wikipedia.org/wiki/Boolean_satisfiability_problem satisfiability]
~~~

~~~
{recent news}
- [https://www.sciencenews.org/article/new-algorithm-cracks-graph-problem graph iso in quasi-poly time ?]
- [https://anthonybonato.com/2017/01/18/graph-isomorphism-and-babais-proof/ Bonato on Babai]
- [https://en.wikipedia.org/wiki/Time_complexity\#Quasi-polynomial_time time complexity]
~~~

~~~
{}{raw}
<h2 id="knapbf">0-1 knapsack: brute force</h2>
~~~

~~~
for each subset
- if weight ok and value best-so-far
-- print subset, update best-so-far-value
~~~

~~~
{}{}
def powerset(iterable): # all subsets
#https://docs.python.org/2/library/itertools.html#recipes
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

def bruteforce(n,val,wt,W):
  indices = range(n) # [0 1 .. n-1]
  bestVal = 0
  for indexset in powerset(indices):
    v, w = 0, 0
    for t in indexset: 
      v+= val[t]
      w+= wt[t]
    if w <= W and v > bestVal:
      bestVal = v
      print indexset, w, v
~~~

~~~
{brute force 0-1 knapsack}
runtime ?
- {{ &Omega;(2<sup>n</sup>) }}

correct ?
- yes: considers every subset

can we do better ? ? ?
~~~

~~~
{}{raw}
<h2 id="knapdp">0-1 knapsack: dynamic programming</h2>
~~~

~~~
subproblems ?
- by set of allowed items
- by allowed total weight

set of allowed items ?
- {{ {1} &nbsp; {1 2} &nbsp;  {1 2 3} &nbsp;  ... }}

allowed total weight ?
- {{ 0 &nbsp; 1 &nbsp; 2 &nbsp; ... }}

{{K[w][j]}}: best value with items {{1..j}}, total wt {{&le;}} w

compute {{K[w][j]}}
- item j out? value {{K[w][j-1]}}
- item j in? value {{K[w-w_j][j-1] + v_j}}
- so set {{K[w][j]}} to max of these two
~~~

~~~
{0-1 knapsack: dp}{}
#python indexing: val[t],wt[t] corresponds to K[][t+1]
# K[w][j]  will be best value knapsack, weight at most w, 
#         using subset of { item_1, item_2, ..., item_j }
# warning: K[][] needs an extra 0th column, so its indices
#   are off-by-one w.r.t. val[], wt[]
# eg. val[v], wt[t]  refer to t+1'st item
#     also K[][t+1]  refers to t+1'st item
def solveknapsack(val,wt,W):
  n = len(val)
  K = [[0 for j in xrange(n+1)] for w in xrange(W+1)]
  for j in range(1,n+1):
    for w in range(W+1):
      K[w][j] = K[w][j-1] if w < wt[j-1] \
        else max(K[w][j-1], K[w-wt[j-1]][j-1]+val[j-1])
~~~

~~~
{0-1 knapsack: dp example}{}
max value knapsack, weight at most 27 ?
 val   [ 6,  9,  7,  9,  8,  7]   
  wt   [11,  6,  8, 10,  8,  9] 
 0 [ 0,  0,  0,  0,  0,  0,  0]
 1 [ 0,  0,  0,  0,  0,  0,  0]
 2 [ 0,  0,  0,  0,  0,  0,  0]
 3 [ 0,  0,  0,  0,  0,  0,  0]
 4 [ 0,  0,  0,  0,  0,  0,  0]
 5 [ 0,  0,  0,  0,  0,  0,  0]
 6 [ 0,  0,  9,  9,  9,  9,  9]
 7 [ 0,  0,  9,  9,  9,  9,  9]
 8 [ 0,  0,  9,  9,  9,  9,  9]
 9 [ 0,  0,  9,  9,  9,  9,  9]
10 [ 0,  0,  9,  9,  9,  9,  9]
11 [ 0,  6,  9,  9,  9,  9,  9]
12 [ 0,  6,  9,  9,  9,  9,  9]
13 [ 0,  6,  9,  9,  9,  9,  9]
14 [ 0   ?   ?   ?   ?   ?   ?]
    
14 [ 0,  6,  9, 16, 16, 17, 17]
15 [ 0,  6,  9, 16, 16, 17, 17]
16 [ 0,  6,  9, 16, 18, 18, 18]
17 [ 0,  6, 15, 16, 18, 18, 18]
18 [ 0,  6, 15, 16, 18, 18, 18]
19 [ 0,  6, 15, 16, 18, 18, 18]
20 [ 0,  6, 15, 16, 18, 18, 18]
21 [ 0,  6, 15, 16, 18, 18, 18]
22 [ 0,  6, 15, 16, 18, 24, 24]
23 [ 0,  6, 15, 16, 18, 24, 24]
24 [ 0,  ?   ?   ?   ?   ?   ?]
 
24 [ 0,  6, 15, 16, 25, 26, 26]
25 [ 0,  6, 15, 22, 25, 26, 26]
26 [ 0,  6, 15, 22, 25, 26, 26]
27 [ 0,  6, 15, 22, 25, 26, 26]
-----------------------------
solution
24 [0, 1, 0, 1, 1, 0] val 26
~~~

~~~
{dp knapsack wc superpoly}
- consider these instances
--  each value and weight has n bits, so in {{[ 2<sup>n-1</sup> ... 2<sup>n</sup> -1 ]}}
-- average weight {{&nbsp; ( 2<sup>n-1</sup>+2<sup>n</sup> -1 ) / 2 &nbsp; &asymp; &nbsp; .75 2<sup>n</sup>}}
-- let W be n times average weight, so {{&nbsp; &asymp; .75 n 2<sup>n</sup>}}
-- input size {{&nbsp; &Theta;( n<sup>2</sup> ) &nbsp; }} bits
- brute force wc time
-- {{2<sup>n</sup>}} subsets, each weight {{&Theta;(n)}} bits
-- for each subset t of size s, perform s add'ns
-- total number of add'ns is sum, over all subsets t, of size(t)-1
-- sum, over all subsets t, of size(t) is {{n*2<sup>n-1</sup>}}
   (pair each subset with its complement)
-- total number add'ns: {{ &nbsp; n*2<sup>n-1</sup> - 2<sup>n</sup> &nbsp; }} in {{&nbsp; &Theta;( n * 2<sup>n</sup>)}}
-- avg add'n involves n/2 weights, sum has {{&nbsp; &Theta;( n<sup>2</sup> )}} bits, time {{&nbsp; &Theta;(n<sup>2</sup>)}}
-- total alg'm time in {{&nbsp; &Theta;( n<sup>3</sup> 2<sup>n</sup> )}}
- dyn. prog. wc time 
-- avg wt. {{&nbsp; &Theta;(n<sup>2</sup>)}} bits
-- avg add'n cost {{&nbsp; &Theta;(n<sup>2</sup>)}}
-- wc W {{&nbsp; &Theta;(n * 2<sup>n</sup>)}}
-- {{&Theta;(n<sup>2</sup> * n * W ) &nbsp; = &nbsp; &Theta;( n<sup>4</sup> 2<sup>n</sup> )}}, worse than BF!
- let {{&nbsp; t = n<sup>2</sup> &nbsp; so n = &radic; t &nbsp; &nbsp; let f(t) = t<sup>2</sup> 2<sup>&radic;t</sup>}}
-- each alg. input {{&nbsp; &Theta;( t ) bits &nbsp; }} wc time {{&nbsp; O( f(t) )}}
- [https://en.wikipedia.org/wiki/Time_complexity\#Table_of_common_time_complexities time complexities]
- {{2<sup>O(log t)</sup>}} polynomial
- {{2<sup>o(t)</sup>}} subexponential 
- {{2<sup>&Omega;(t)</sup>}} exponential 
- here
-- {{ lim<sub>t &rarr; &infin;</sub> &nbsp; 2<sup>k log t</sup> / f(t) &nbsp; = 0, &nbsp;}} so f(t) superpoly
-- {{f(t)/2<sup>t</sup> &nbsp; = &nbsp; t<sup>2</sup> 2<sup>&radic;t</sup> &nbsp;/&nbsp; 2<sup>t</sup> &nbsp; = &nbsp; t<sup>2</sup> (2 / 2<sup>&radic;t</sup>)<sup>&radic;t</sup>}}
-- {{ lim<sub>t &rarr; &infin;</sub> &nbsp; f(t) / 2<sup>t</sup> &nbsp; = 0, &nbsp;}} so f(t) subexp
~~~

~~~
{}{raw}
<h2 id="knapdpv">0-1 knapsack: another DP solution</h2>
~~~

~~~
- how do you like DPW if
-- each weight is in {{ { 2^{n-1}, ..., 2^n -1 } }}  and
-- each values is in {{ { 1, ..., 2n } }} ?
- organize subproblems by value (instead of weight)
- also useful in constructing knapsack PTAS (polytime approx scheme)
~~~

~~~
{}{}
def knapDPV(val,wt,V): #dynamic programming by value
  # A[v][j] is min weight of subset of items 0..j with exact val v
  infinity = 1 + sum(val) #larger than any possible sum of values
  n = len(val)
  A = [[infinity for j in xrange(n)] for v in xrange(V+1)]
  for j in range(n): A[0][j] = 0
  A[val[0]][0] = wt[0]  # end initialization
  for v in range(1,V+1):  # row    0 already initialized
    for j in range(1,n):  # column 0 "       "
      A[v][j] = A[v][j-1] if v < val[j] \
        else min(A[v][j-1], A[v - val[j]][j-1] + wt[j])
~~~

~~~
- [http://webdocs.cs.ualberta.ca/~hayward/304/asn/knapdpvdemo.pdf
  knapsack DPV hand-trace]
~~~

~~~
{knapsack DPV}
- time in O(nV), but V can be big
- correctness ?
~~~ 

~~~
{}{raw}
<h2 id="setcover">setcover</h2>
~~~

~~~
{brute force set cover}{}
from random    import randint
from itertools import chain, combinations

def powerset(iterable):
#https://docs.python.org/2/library/itertools.html#recipes
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

def settostring(s,n): # set as pretty binary string
  x = ''
  for j in range(n):
    if j in s:   x += ' *'
    else:        x += ' -'
  return x

def gensubset(a,b,n): # subset of 0..n-1, each with prob a/b
  s = set()
  for j in range(n): 
    if randint(1,b) <= a: s.add(j)
  return s

def initsubsets(n,m,a,b): # m subs of 0..n-1
  # an element in a subset with probability a/b
  L = [] # list of subsets
  for _ in range(m): L.append(gensubset(a,b,n))
  sumlen = 0
  for x in L:
    sumlen += len(x)
    print settostring(x,n)
  print 'universe size', n,' avg subset size', sumlen*1.0/m
  return L

def bruteforce(n,m,L):
  indices = range(m) # [0 1 .. m-1]
  bestcover = ()
  for indexset in powerset(indices):
    subsetunion = set()
    for t in indexset: 
      subsetunion |= L[t]
    # indexsets sorted by cardinality, so first found is min
    if len(subsetunion) == n: 
      print 'size', len(indexset), 'cover', indexset
      for j in indexset: print settostring(L[j],n)
      return # so return once one is found

n,m,a,b = 40,25,2,11
L = initsubsets(n,m,a,b)
bruteforce(n,m,L)
~~~

~~~
{set cover example}{}
set     universe
    0 1 2 3 4 5 6 7 8 9
0   - - * - - - - * * -
1   - * - - - - - * - *
2   - - - - * - - * * -
3   - - - - - - - - * -
4   - - - - - * - * * -
5   - - - - * - * * - *
6   * * - * - - - - - -
7   - - - - - - * * - -

min cover (0, 4, 5, 6)
0   - - * - - - - * * -
4   - - - - - * - * * -
5   - - - - * - * * - *
6   * * - * - - - - - -
~~~

~~~
{another example}{}
set                    universe
    0 1 2 3 4 5 6 7 8 910111213141516171819
0   * - * - - - * - * - - - - * - * - - * -
1   - - - * - * - * - - - * - - * * - - - -
2   - * - - - - * - * - - - - - - - - * - -
3   - - - * - - - - - - - - - - - - - - * -
4   * - - - - * - * - * - - - - - * - - - -
5   - - - * - - * * - - - - - * - * - - - -
6   - - - - - * - - * - - * * - - * - - * -
7   - * - - * - - - * - * * * - * * - - - -
8   - * - - * - - - * - - - - * * - - - * *
9   * - * - - - * * - - - - - - - - - - - *
10  - * - - - - * - - * * - * - - - * - * -

min cover (0, 1, 2, 8, 10)
0  * - * - - - * - * - - - - * - * - - * -
1  - - - * - * - * - - - * - - * * - - - -
2  - * - - - - * - * - - - - - - - - * - -
8  - * - - * - - - * - - - - * * - - - * *
10 - * - - - - * - - * * - * - - - * - * -
~~~

~~~
{greedy set cover}
- [https://webdocs.cs.ualberta.ca/~hayward/304/asn/greedysc.pdf 
   greedy set cover]
- Dasgupta et al. section 5.4
~~~

~~~
{school location problem}{}
  - input: some towns, and a threshold distance (e.g. 30 km)
  - problem: build schools in minimum number of towns,
      s.t. distance to school at most the threshold
  - model as a graph problem G = (V,E)
      - V = set of towns
      - E = pairs of towns at most threshold distance apart
      - for each town, its neighbor-set is set of all towns
          whose distance away is at most the threshold  
      - problem: pick the minimum number of neighbor-sets that cover V
      - in other words: this is an instance of min set cover
~~~


~~~
{how hard is set cover?}{}
- answer: NP-hard, by reduction from vertex cover
~~~
     
~~~
{}{raw}
<h2 id="sat">satisfiability</h2>
~~~

~~~
given a boolean formula in CNF, is it satisfiable?

assume input has n variables and m clauses

brute force: {{2<sup>n</sup>}} possible assignments

brute force: worst case time {{ O(m2<sup>n</sup>) }}
~~~

~~~
{sat example}{}
( x1  or  ~x2  or  x3 )  and
( x1  or  ~x2  or  x4 )  and
( x1  or   x2  or  x3 )  and
( x1  or   x2  or  x4 )  and
(~x1  or  ~x2  or ~x3 )  and
(~x1  or  ~x2  or  x3 )  and
(~x1  or   x2  or  x4 )  and
( x2  or  ~x3  or ~x4 )  and
( x2  or   x3  or ~x4 )

[  1, -2,  3]
[  1, -2, -4]
[  1,  2,  3]
[  1,  2,  4]
[ -1, -2, -3]
[ -1, -2,  3]
[ -1,  2,  4]
[  2, -3, -4]
[  2,  3, -4]

satisfiable (0, 1, 1, 0)
~~~

~~~
{another example}{}
[ -4, -5,  6]
[ -4,  5, -6]
[ -2,  4, -5]
[ -2,  5, -6]
[ -1, -3, -4]
[ -1, -3,  4]
[ -1,  4, -6]
[  1, -4,  5]
[  1, -3, -5]
[  1, -3,  4]
[  1,  5, -6]
[  2,  5,  6]
[  3, -5, -6]
[  4, -5,  6]
[  4,  5, -6]
[  4,  5,  6]

satisfiable (1, 1, 0, 1, 0, 0)
~~~

~~~
{brute force sat solver}{}
# I/O: literals x1 -x1 x2 -x2 ... rept'd by 1 -1 2 -2 ...
def clause_sat(clause,asnmt):
  for var in clause:
    if ((var>0 and asnmt[  var-1 ] == 1) or
        (var<0 and asnmt[-(var+1)] == 0)):
      return True
  return False

def sat(formula, asnmt):
  for clause in formula:
    if not clause_sat(clause,asnmt): return False
  return True

n, k, m = 20, 3, 90  #max m: (n choose k)(2^k)
myf = formula(n,k,m)
show(myf, False)

allAsn = product([0,1],repeat=n)
for a in allAsn:
  if sat(myf,a): print 'satisfiable', a
~~~