Summary Algorithms and Data Structures (c) 2004-02-24 Fabian M. Suchanek http://www.mpi-inf.mpg.de/~suchanek/personal/texts/summaries/ads.txt This is a summary of the course "Algorithms and Data Structures" held by Priv.-Doz. Dr. Peter Sanders and Dr. Kavitha Telikapalli in the WS 2003 at Saarland University. Please note that, although I am also a tutor of this course, this is a student's summary and not official. I would like to thank the participants of my tutorial, who helped finding the bugs in this summary, especially David Steurer, Gernot Gebhard and Peter Kolloch. Furthermore, I owe thank to Thomas Schultz, who scanned the largest part of the summary rigorously for errors. ___ | Since this text contains some ASCII-art, it should only / \ | / be viewed or printed with a fixed font (eg Courier New). ( ) |< If the font is all right, an "OK" appears here: \___/ | \ By reading the following text, you accept that the author does not accept any responsibility for the correctness or completeness of this text. If you have any corrections or remarks, please send me a mail. This is the only way to make the publication of this summary useful for me, too. My e-mail address is f.m.suchanek@zweb.de, but the letter 'z' has to be removed from the address. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Prerequisites ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Set: An unordered collection of different elements. S = {s1,...sn} S = { x | p(x) } is the set of all those x for which p(x) holds Epsilon of a set S: An element of S. Notation: eps(S) := x, such that x in S eps({x}) := x // This operator has been invented by Hilbert. Since it is // the only way to convert a set to an element, this summary makes // use of it. Size, cardinality of a set: The number of its elements. Notation: |{s1,s2,...sn}| = n Intersection of two sets: The set of those elements which occur in both of them. Notation: X /\ Y Union of two sets: The set of all of their elements. Notation: X \/ Y Difference of two sets: The set of those elements which occur in the first one, but not in the second one. Notation: X \ Y Disjoint sets: Sets, the intersection of which is empty. Picking: Taking an element of a set and deleting it from the set pickAny(S) := (z:=eps(S), S:=S\{z}, z) Multiset: An unordered collection of elements. Unlike a set, a multiset may contain the same element more than once. Tupel: An ordered collection of elements. A tuple can be assigned to a set, but not vice versa. Notation: t=(t[1],t[2],...t[n]) Example: t=(1,'A',abc) is a tupel Dimension of a tupel: The number of its elements. Pair: A tupel of dimension 2. Triple: A tuple of dimension 3. Subset of a set S: A set, the elements of which are all contained in S. Notation: S' =< S Superset of a set S: A set, which contains all elements of S. Notation: S' >= S Partition of a set S: A set of subsets S[1],...S[n] of S, such that each element of S is in exactly one of them. Cartesian Product of two sets A and B: The set of all possible pairs, the first element of which is in A and the second element of which is in B. Notation: A x B = C Example: {1,2} x {a,b,c} = { (1,a), (2,a), (1,b), (2,b), (1,c), (2,c) } n-th Power of a set S: The set S^n := S x S x S x ... S where 'S' occurs n times. Relation over two sets A and B: A subset of the cartesian product of A and B. Notation: * (a,b) is an element of the relation R * r(a,b) * a r b * (a,b) is in R Example: The relation ">" is the set of all pairs of real numbers, in which the first number is greater than the second. It holds: ">" = { (1,0), (27.3,12.5), (17,-12), ... } One writes: >(3,2) or 3 > 2 or (3,2) is in ">". Reflexive relation: A relation R < A x A such that R(a,a) for all a in A. Symmetric relation: A relation R < A x A such that if R(a,b), then R(b,a) for all a,b in A. Antisymmetric relation: A relation R < A x A such that if R(a,b) and R(b,a) then a=b for all a,b in A. Asymmetric relation: A relation R < A x A such that if R(a,b), then not R(b,a) for all a,b in A. Complete relation: A relation R < A x A such that R(a,b) or R(b,a) for all a,b in A. Transitve relation: A relation R < A x A such that if R(a,b) and R(b,c), then R(a,c) for all a,b,c in A. Strict relation: A relation R < A x A such that not R(a,a) for all a in A. Linear order on a set A: A relation R < A x A, which is asymmetric, transitive and complete. Usually, this relation is denoted by '<'. Example: '<' is a linear order on the set of natural numbers. Simple order on a set A: A relation R < A x A, which is reflexive, antisymmetric, transitive and complete. Usually, this relation is denoted by '=<'. Example: '=<' is a simple order on the set of natural numbers. Domain of a relation over A and B: A. Range of a relation over A and B: B. Function: A relation f over two sets A and B, such that for every a in A there is either none or exactly one element b in B with a f b. Notation for f: f:A->B Notation for "a f b": f(a)=b. Notation for the set of all functions f:A->B: B^A There are |B^A| = |B|^|A| many of these functions f:A->B. Result of a function f:A->B for an argument x: The element of B such that f(x)=b. Real: The set of real numbers. For a detailed definition and associated operations, cf "Constraints over the Reals" (cr.txt) Real+: The set of positive real numbers with 0. Integer: The set of integer number. For a detailed definition and associated operations, cf "Constraints over the Reals" (cr.txt) Natural: The set of positive integers, including 0. Interval with the real limits a and b: The set [a,b] := {x | a =< x =< b}. Prime: A natural number greater than 1 which can only be divided by 1 and itself. Order of a function f:Natural->Real: A function g:Natural->Real, such that there is x0 in Real and c in Real, such that f(x)x0. Class of an Order f:Natural->Real: The set of functions having f as an Order. Notation: O(f) Notation for "g:Natural->Real is an element of O(f)": g in O(f) Necessity of a state A for a state B: The fact that B cannot be without A, B => A. Sufficiency of a state A for a state B: The fact that, whenever A is, B is also, A => B. In order to prove the equivalence of two states, it is sufficient to prove necessity and sufficiency. iff: if and only if. Operation, Algorithm: A sequence of instructions, which may have a result. Recursive algorithm: An algorithm X, one instruction of which is the execution of X itself. Termination of an algorithm: Its property that it cannot be executed forever. Partial correctness of an algorithm: Its property that it gives the right result, if it terminates. Partial correctness does not yet guarantee that the algorithm terminates. Total correctness of an algorithm: Its property that it is partially correct and terminates. In this case, the algorithm does really what it is intended to do. Time of an algorithm for an input size n: The order of the function which calculates the time needed to execute the algorithm, in dependence upon the input size n. One says "The algorithm is in time O(g)". Polynomial in x: An expression of the form a[0]*x^0 + a[1]*x^1 + ... a[n]*x^n where x is a real variable and the a[i] are real constants. // See "Constraints over the Reals" (cr.txt) for more details Algorithm with a polynomial time: An algorithm, the time of which can be expressed by a polynomial in the input size. NP-hard problem: A problem, for which no algorithm is known, which solves the problem in polynomial time. Space of an algorithm for an input size n: The order of the function which calculates the number of cells needed to execute the algorithm, in dependence upon the input size n. One says "The algorithm is in space O(g)". RAM, Random Access Machine: A model of a computer with an infinite amount of cells. A RAM can execute certain instructions. For a detailed definition cf "Theoretische Informatik" (infod.txt, German) Primitive type: A set of finite size. Boolean: The primitive type {true, false}. Class, Type: A primitive type or a tuple of types together with a tuple of names. Usually, there are operations associated to types. Notation for tuple types: ( : , : , ... ) or CLASS EXTENDS : ... PROCEDURE (...) ... ... The operations can only be called with an instance (s.b.). This instance is referred to as in the procedure. For convencience, the string "this." can be left away when accessing components of . The EXTENDS-clause says that the type inherits all components and operations of the indicated type. If a class inherits an operation and also defines it, then the newly defined operation is used for all instances of the extending class. Data structure, Instance of a type t: If t is a primitive type, then an element of t, if t is a tuple of types, then a tuple of data structures of the elements of t. Notation for the component named c of a datastructure x: x.c Notation for the operation named f of a datastructure x: x.f Notation for a datastructure of type t: (:, :, ...) or (, , ...) Notation for a datastructure of type t constructed from a datatstructure x of some other type: t(x) 1-dimensional array, array, sequence of type t: A type, which is a tuple of t's. Notation for data structures: s= Associated operations: * s[i] := * s.length() := * s1 + s2 := * s[i..j] := nil: A value, which is compatible with all types. Singly linked list of type t: The type (car:t, cdr:list of type t). That is: A list data structure is a pair of a data structure and a list data structure. CLASS List of T car : T cdr : List of T FUNCTION length() : Integer RETURN (this==nil) ? 0 : 1+this.cdr.length PROCEDURE insert(i : Integer, x : T) IF (i==0) THEN this:=List(x,this) ELSE this.cdr.insert(i-1) PROCEDURE pushFront(x : T) this.insert(0,x) PROCEDURE pushBack(x : T) this.insert(this.length,x : T) FUNCTION popFront() : T RETURN this.delete(0) FUNCTION popBack() : T RETURN this.delete(this.length-1) FUNCTION last() : T RETURN (this.cdr==nil) ? this.car : this.cdr.last() FUNCTION delete(i : Integer) : T IF i==0 THEN x:=this.car; this:=this.cdr; RETURN x ELSE RETURN this.cdr.delete(i-1) FUNCTION delete(x : T) : T IF this==nil RETURN nil IF this.car==x THEN RETURN this.delete(0) RETURN this.cdr.delete(x) FUNCTION get(i : Integer) : T IF i==0 THEN RETURN this.car RETURN this.cdr.get(i-1) PROCEDURE +(l : List of T) this.cdr:=l List of type t: A type, which has operations equivalent to those of a singly linked list of type t. Doubly linked list of type t: The list CLASS doublyLinkedList of type T element : T prev : doublyLinkedList of type T succ : doublyLinkedList of type T with the appropriate operations. List-Stack of type t: The type CLASS ListStack of T EXTENDS List of T FUNCTION pop() : T RETURN this.delete(0) FUNCTION push(x : T) RETURN this.insert(0,x) FUNCTION top() : T RETURN this.car FUNCTIOn empty() : Boolean RETURN this.length==0 Stack of type t: A type, which has operations equivalent to those of a list-stack of type t. Naive Queue of type t: A list of type t with the following operations CLASS NaiveQueue of T EXTENDS List of T this.insert(this.length,x) FUNCTION dequeue(), q.popFront() : T RETURN this.delete(0) FUNCTION head() : T RETURN this.car FUNCTION empty() : Boolean RETURN this.length==0 Queue of type t: A type, which has operations equivalent to those of a naive queue of type t. Naive Priority Queue of type t and priority f:t->Real: The type CLASS NaivePQueue of T with f:t->Real EXTENDS List of T PROCEDURE insert(x) IF this==nil || f(this.car)>f(x)) THEN this:=NaivePQueue(x,this) ELSE this.cdr.insert(x) FUNCTION deleteMin() : T RETURN this.delete(0) PROCEDURE updateKey(x : T), decreaseKey(x : T) this.delete(x) this.insert(x) FUNCTION top() : T RETURN this.car Priority Queue of type t and priority f:t->Real: A type, which has operations equivalent to those of a Naive Priority Queue of type t and priority f. Naive heap of type t and priority f:t->Real: A naive priority queue, which supports the following additional operations: CLASS NaiveHeap of T with f:t->Real EXTENDS NaivePQueue of T FUNCTION pred(x : T) : T IF this==nil THEN RETURN nil RETURN this.cdr.car==x?this.car:this.cdr.pred(x) FUNCTION succ(x : T) : T IF this==nil THEN RETURN nil RETURN this.car==x?this.cdr.car:this.cdr.succ(x) PROCEDURE swap(x,y : T) // only neighbors IF this==nil THEN RETURN IF this.car==x && this.cdr.car==y THEN (l.car,l.cdr.car) := (l.cdr.car, l.car) ELSE this.cdr.swap(x,y) FUNCTION delete(x : T) : T IF (this==nil) THEN RETURN nil IF this.car==x THEN this:=this.cdr RETURN x ELSE RETURN this.cdr.delete(x) FUNCTION locate(x : T) : T // Returns x or next smaller element IF this==nil RETURN nil IF this.car==x || f(this.cdr.car)>f(x) THEN RETURN this.car RETURN this.cdr.locate(x) Heap of type t and priority f:t->Real: A type, which has operations equivalent to those of a Naive Heap of type t and priority f. Reverse of an array X: The array reverse(X) := . Matrix of type t: A sequence of n sequences of type t. The inner sequences all have to be of the same length m. The matrix is then called an n x m matrix. A vector can be seen as n x 1 matrix. Transposed matrix of a n x m M: The n x m matrix t(M) with t(M)[j][i] = M[i][j] Determinant of a 3 x 3 matrix A of type Real: The value det(A) := a*e*i + b*f*g + c*d*h - c*e*g - b*d*i - a*f*h if a b c A = d e f g h i That is: Sum up the products of all \-Diagonals, subtract the products of all /-Diagonals. // For details, see "Constraints over the Reals" (cr.txt) Position, Index of an element in an array data structure: The number of its predecessors in the array. Vector: The type of a tuple of reals. Its operations are defined in the following. // For a detailed definition, see "Algebra" (algebra.txt, German) // or "Analyse de Donnees" (maths.txt, in French) Sum of two vectors x and y: The vector x+y := (x[0]+y[0], ..., x[x.length-1]+y[x.length-1]). r-Multiple of a vector x: The vector vector r*x := (r*x[0], ..., r*x[x.length-1]). Negative of a vector x: The vector -x := (-x[0],...-x[x.length-1]). Absolute value of a vector: The root of the sum of the squares of its components. Notation: |x| := sqrt(SUM {x[i]^2 | i=1..x.length-}) Dot-product of two vectors x and y: The real value x * y := SUM { x[i]*y[i] | i=1..x.length-1 } The dot-product measures the "similarity" of the two vectors. It is 0 if the vectors are orthogonal. // For details, cf "Analyse de Donnees" (maths.txt, in French) n-dimensional array of type t: An array, the elements of which are (n-1)-dimensional arrays of type t. Notation for a[i[0]][i[1]]...[i[n]]: a[i] for an integer vector i Binary representation of a natural number x: The queue binrep(x) := (x==0) ? <0> : binrep(x div 2).enqueue(x mod 2) Natural number of a binary representation b: The natural number unbin(b) := _unbin(b,1) _unbin(b,i) := (b==nil) ? 0 : i*(b.last() mod 2)+unbin( (b.delete(b.length-1),b), i*2 ) Bit: An element of {0,1}. xor: A function mapping two natural numbers a and b to the natural number of the binary representation, which is achieved by setting those bits to 1, for which the binary representations of a and b differ. It holds: a xor a = 0 a xor b = c => a = b xor c a xor b = b xor a a xor a xor b = b Random variable: A function, which maps individuals to real numbers. A random variable represents a characteristic of the individuals. For each individual, the random variable "knows" the degree of markedness of this characteristic. Example: Consider the age of a person. A random variable might be A, such that A(grand_ma)=99, A(fabian)=23 etc // For a more detailed definition see "Statistik" (statistik.txt, // German) Expectation of random variable X: The sum of the results of X for all possible individuals, divided by the total number of all possible individuals. This is a kind of theoretical mean value. Notation: E(X)=y Linearity of expectation: The fact that E(X + Y) = E(X) + E(Y). Probability of a random variable X taking the value y: The number of all possible individuals i, for which X(i)=y, divided by the total number of all possible individuals. Notation: P(X=y) Probability of X=x if Y=y: The probability that the random variable X takes the value x, provided that the random variable Y takes the value y. Notation: P(X=x | Y=y) Independant random variables: Two random variables X and Y, such that P(X=x | Y=y) = P(X=x). That is: Whether Y=y or not, this does not influence P(X=x). Probability of X=x and Y=y: The probability that the random variable Z(i) := (X(i)==x and Y(i)==y) ? 1 : 0 takes the value 1. It holds P(X=x and Y=y) := P(X=x)*P(X=x | Y=y) Probability of X=x or Y=y: The probability that the random variable Z(i) := (X(i)==x or Y(i)==y) ? 1 : 0 takes the value 1. It holds P(X=x or Y=y) := P(X=x)+P(Y=y) - P(X=x and Y=y) k-th harmonic number: The sum H(k) := 1/1 + 1/2 + 1/3 + ... + 1/k It holds that H(k) < ln(n)+1. Set sum: The function SUM, which maps a set of real numbers to their sum. Notation: SUM {x[1],x[2],...x[n]} := x[1]+x[2]+...+x[n] Lemma, Theorem: A fact. Sorted sequence of a set A: A sequence S of type A, which contains all elements of S, such that S[i] Q // Transition function FUNCTION runOn(S : array of I) : boolean state := s0 pos := 0 WHILE pos<|S| DO state := d(state,S[pos]) pos := pos + 1 RETURN (state in F) A FSA "reads" a string symbol by symbol. At each symbol read, the state is changed according to the function d. If the string is over and the FSA remains in a state marked as "final", then the FSA is said to "accept" the string. A FSA may be depicted as a graph where every node stands for a state and the edges show a transition from one state to another, given a particular input symbol. A FSA may also be depicted by a transition table, which shows the new state(s) for a particular current state and a particular input. // For details, see "Computational Linguistics" (cl.txt) or // "Theoretische Informatik" (infod.txt, in German) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Graphs ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ See also "Graphes et Complexite" (gc.txt, in French) Graph: The type CLASS graph vertices : Set // often noted V edges : Subset of vertices x vertices // often noted E A graph is usually seen as consisting of "nodes" (the elements of V) and "connections" between them (the elements of E). Graphs are usually drawn as networks, in which the nodes are connected by arrows. In the following, V is always the set of nodes and E is always the set of edges. Graphs usually have the following operations: * Access to an edge * Find outgoing edges of a node * Find an edge * Insert edge * Delete edge * Insert node * Delete node Vertex, Node of a graph (V,E): An element of V. The number of nodes is denoted by n. Usually, the nodes are simply integers 0..n-1. Edge of a graph (V,E): An element of E. The number of edges is denoted by m. For an edge e=(u,v), one says that e goes from u to v. Sometimes the edge (u,v) is noted by "u->v". Usually, the edges are numbered. We define CLASS edge startnode : Node endnode : Node Reverse edge of an edge (u,v): The edge (v,u). Going to a node, looking at a node: Shifting attention towards this node. Traversing an edge, following an edge: Going from its start-node to its end-node. Successor of a node u: A node v, for which there is an edge (u,v). Outdegree of a node v: The number of edges going out of v. outdegree(v) := |{(v,u) in E}| Indegree of a node v: The number of edges going to v. indegree(v) := |{(u,v) in E}| Self-loop: An edge, the start-node of which is its end-node. Usually, self-loops are disallowed. Bidirected graph: A graph, in which for every edge (u,v), there exists one edge (v,u). Undirected graph: A bidirected graph, in which the edges (u,v) and (v,u) are considered equal. This is why the edge is usually noted by {u,v}. An undirected graph is usually drawn with simple lines instead of arrows between the nodes. Degree of a node v: The number of edges involving v. If the graph is directed, then it holds degree(v) := indegree(v) + outdegree(v) In a graph (V,E), the number of nodes with an odd degree is even: SUM {degree(s) | degree(s) mod 2==1} + SUM {degree(s) | degree(s) mod 2==0} = SUM {degree(s)} = 2*|E| 2*|E| is even SUM {degree(s) | degree(s) mod 2==0} is even => SUM {degree(s) | degree(s) mod 2==1} is even => |{degree(s) | degree(s) mod 2==1}| is even Sub-graph of a graph (V,E): A graph (V',E') with V'=Real. These edge costs can be thought of as an amount of money to be paid when the edge is traversed. Notation: c(u,v) := c((u,v)) for an edge (u,v) Weighted graph: A graph with edge weights. We define CLASS WeightedGraph EXTENDS Graph c:FUNCTION edges->Real Path in a graph (V,E): A sequence of nodes , such that v[i] in V and (v[i],v[i+1]) in E. A path can also be given as the sequence of the respective edges: <(v[0],v[1]), (v[1],v[2]), ... (v[k-1],v[k])>. Usually, a path is noted by v[0]->v[1]->...v[k] or as v[0]~~>v[k]. Start node of a path: Its first element. End node of a path: Its last element. Node reachable from a node s: A node t, for which there is a path from s to t. Length of a path in a graph: The number of its nodes minus 1. Notation: length() := k Weight, Cost of a path p in a weighted graph: The sum of the costs of the edges of p. Notation: c() := SUM {c(v[i],v[i+1]) | i=0..k-1} Simple path: A path, the nodes of which are all distinct. Walk: A path, which is not simple. Bottleneck weight of a path p in a weighted graph: The minimum of the weights of the edges of p. Cycle: A path, the start node of which is at the same time its end node. Simple cycle: A cycle, the nodes of which are all different except for the start and end node. Hamiltonian cycle: A simple cycle which contains all nodes. Finding a Hamiltonian cycle is an NP-complete problem. Euler cycle: A cycle which contains all edges once. Connected graph: A graph, in which there is for two nodes always a path from one to the other, if one adds to each edge (u,v) its reverse edge (v,u). This definition simply says that the graph does not consist of separated components. Strongly connected graph: A connected graph, in which there is a path from every node to every other node. Complete graph: A graph (V,V x V). That is: Every node is connected to every other node. Eulerian graph: An undirected graph which contains an Euler cycle. A graph is eulerian, iff * it is connected and * every node has even degree Necessity: * If the graph was not connected, there could not be a cycle touching all its edges. Therefore, the graph must be connected. * If the Euler cycle enters into a node, it must also quit this node again. Since every edge is only to be visited once, each passing of a node requires 2 separate edges. Sufficiency: By the following algorithm, an Euler cycle can be found if the above conditions are met: FUNCTION eulerCycle_(VAR (V,E):Graph,s: V): List of E IF s==nil THEN s := eps(V) result : List of E := <> u : V := s REPEAT next : V := eps({v | (u,v) i E}) List.pushBack((u,next)) E := E \ {(u,next)} u := next UNTIL u==s WHILE Ex i: degree(result.get(i).endNode) mod 2==1 DO s := result.get(i) small : List of E := eulerCycle((V,E),s.endNode) small.cdr := s.cdr; s.cdr := small RETURN result FUNCTION eulerCycle((V,E):Graph): List of E RETURN eulerCycle_((V,E):Graph,eps(V)) Start at any node. Follow an edge to another node and delete this edge. Continue this way until you get to a node with no more edges. Since every node has an even degree, each node can be entered and left. If one gets to a node with no exit, then this must be the start node, because it is the only one with an uneven degree. Start at another node and do the same. Continue until no more edges are left. This way, many cycles are created. Every one of them has at least one node in common with another cycle. If this were not the case, the graph would not be connected. These common nodes can be used to glue together the cycles, creating one large cycle -- the Euler cycle. The algorithm runs in time O(|E|). acyclic graph: A graph which does not contain any cycles. Undirected tree: An undirected graph, in which one of the following properties holds: 1. There is exactly one path from one node to another 2. The graph is connected and has exactly n-1 edges 3. The graph is connected and acyclic. These properties are equivalent: * 1 => 2 If there is a path from one node to another, then the graph is connected. If there is a path from one node to another, then the graph cannot have less than n-1 edges, since every node needs one edge. If the graph had more than n-1 edges, it would have a superfluous edge, which could be used to construct an alternative path. * 2 => 3 If the graph is connected, then it is connected. The graph has n-1 edges. Suppose it contains a cycle. Then we can delete one edge without destroying the connectedness. Then the graph contains n-2 edges, which can never connect all nodes. This contradiction implies that the graph cannot contain a cycle. * 3 => 1 Since the graph is connected, there exists one path from one node to another. Suppose there were two paths from one node u to another node v. Then u~~>v~~>u would be a cycle. Since 3 says that the graph does not contain any cycle, each path must be unique. Forest: An undirected acyclic graph. A forest can be seen as consisting of multiple trees (which is quite a reasonable definition). Such a graph can at most have n-1 edges. If it had more, it would contain a cycle, see definition of "undirected tree". Directed tree: A directed graph, which has a node r, such that there is exactly one path from r to every other node. This entails * Connectedness If the graph was not connected, there could not be this path from r to a node in another component. * No cycles If the graph allows a cycle around a node u, then there are infinitley many paths from r to u, namely r~~>u, r~~>u~~cycle~~>u, etc. * At most one incoming edge Every additional incoming edge would allow another path from r. If all edges are reversed in such a graph, then this is also called a tree. Root of a directed tree: The node, from which there is a path to every other node. Parent of a node u in a directed tree: The node for which there is an edge to u. The parent of the root is the root. Child of a node u in a directed tree: A node v for which there is an edge (u,v). A node may have multiple children. Leaf of a directed tree: A node which has no children. Interior node of a directed tree: A node which has children. Ordered directed tree: A directed tree, in which each node maintains an ordering of its children. DAG, directed acyclic graph: A graph without cycles. This type of graphs is for instance used in multiple inheritance "trees". A graph can be tested for being a DAG by the following algorithm: FUNCTION isDAG((V,E) : Graph) : Boolean WHILE Ex v in V: outdegree(v)==0 DO V := V \ {v} E := E \ {(u,v} | u in V} \ {(v,u} | u in V} RETURN (V=={}) That is: While there is still a node without outgoing edges, delete this node. Iff no node remains, then the graph is a DAG. This is because if the graph contains a cycle, then this cycle is supported by outgoing edges. Each node protects its predecessor from being deleted. // For an application of DAG, see "Representing Ontological Structures // in CLOS, Java and FAST" (ba.doc) Graph representation: A type, which provides the operations of graph in an efficient way on a computer. // Since the operations can mostly be implemented trivially, // this summary only gives their run-times Edge sequence representation: The graph representation CLASS edgeSequence edges : List of Edge nodes : List of Node This class uses a list of all edges <(u1,u2),(u5,u2), ...> and a list of all nodes of the graph. (The node array is necessary to store nodes which have no edges.) Operations with their time: * Access to an edge: O(m) * Find outgoing edges of a node: O(m) * Find an edge: O(m) * Insert edge: O(1) * Delete edge: O(m) (first find the edge) * Insert node: O(1) * Delete node: O(n) Additional information for edges or nodes can be stored in partner arrays of the two arrays. Space: O(m+n) Adjacency array: The graph representation CLASS adjArray edgeArray : array of Node nodeArray : array of Integer This class uses two arrays: The "edge array" stores, for each node u, the end nodes of all edges going out of u. The "node array" stores, for each node u, the index of the cell in the edge array, which contains the first outgoing edge of u. The last element of the node array stores the length of the edge array. Operations with their time: * Find outgoing edges of a node: O(1) * Find an edge: O(m) * Insert edge: O(m+n), needs to create a new edge array * Delete edge: O(m), a second node array can be used, which stores the index of the last edge of a node. * Insert node: O(n), needs to create a new array * Delete node: O(n), the node array needs to be re-arranged Additional information for edges and nodes can be stored in partner arrays of the two arrays. Space: O(m+n) Adjacency list: The graph representation CLASS adjList list : array of List of Node Here, each array element with index i stores a list of the successors of node i. Operations with their time: * Find outgoing edges of a node: O(1) * Find an edge: O(m) * Insert edge: O(1) * Delete edge: O(m) (first find the edge) * Insert node: O(n), needs to create a new array * Delete node: O(n), the array needs to be re-arranged Space: O(m+n), but three times as much as adjacency arrays. Adjacency matrix: The graph representation CLASS adjMatrix matrix : array of array of Bit The matrix contains a 1 at position i,j iff there is an edge from node i to node j. Operations with their time: * Access to an edge: O(1) * Find outgoing edges of a node: O(n) * Find an edge: O(1) * Insert edge: O(1) * Delete edge: O(1) * Insert node: O(n^2), needs to create a new array * Delete node: O(n^2), needs to create a new array Space: O(n^2) If an adjacency matrix is raised to the k-th power, the position i,j indicates the number of paths of length k between nodes i and j. For a definition of k-th power and a proof of the path property, see "Graphes et Complexite" (gc.txt, French). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Graph traversal ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Graph traversal: The process of systematically going through all edges and nodes of a graph. This process starts in a certain start node. Shortest path from node u to node v: The path with the start node u and the end node v, which has the smallest length of all possible such paths. Shortest sub-path lemma: The fact that if is the shortest path from u to v, then is the shortest path from u to v'. Proof: Suppose there was a shorter path from u to v'. Then this path could be used to go from u via v' to v. This path would be shorter than the original one, leading to a contradiction. BFS, Breadth-first search: A graph traversal, which first looks at the start node s, then at all nodes reachable from s by one edge, then all nodes reachable from s by two edges etc. Algorithm: PROCEDURE bfs((V,E) : Graph, s : Node) q : Queue of Node := marked : Array[|V|] of Boolean := marked[s] := true WHILE q != nil DO u := q.popFront() FOREACH (u,v) in E DO IF marked[v]==false THEN q.pushBack(v) marked[v]:=true The idea of this algorithm is to maintain a queue, which contains the nodes to be explored next. Since the children of an explored node are enqueued at the end of the queue, it is guaranteed that the nodes of one level are finished before starting with the next level. This algorithm only works for nodes which are reachable from s. The algorithm runs in time O(m+n), because every edge and every node is explored exactly once. DFSAction: The following type: CLASS DFSAction PROCEDURE traverse((u,v) : Edge) // To be executed if the edge (u,v) is traversed PROCEDURE root(s : Node) // To be executed if the node s demarks a new start node PROCEDURE backtrack((u,v) : Edge) // To be executed if the edge (u,v) is traversed in the // opposite sense when there is nothing interesting at node v Such a DFSAction can be plugged into a graph traversal algorithm. The algorithm will execute the operations of the DFSAction in the order it traverses the graph. // This is the object-oriented interpretation of the blue lines // in the algorithms in the slides DFS, Depth-first search: A graph traversal, which first looks at the start node and then does depth-first searches for all children of the start node. Algorithm: PROCEDURE dfs((V,E) : Graph, s : Node, a : DFSAction) S : Stack of Node := nil markedNode : Array[|V|] of Boolean := markedEdge : Array[|E|] of Boolean := FOREACH s in V DO IF markedNode[s]==false THEN markedNode[s] := true S.push(s) a.root(s) WHILE S != nil DO v := S.top() IF Ex e=(v,w) in E : markedEdge[e]==false THEN a.traverse(v,w) markedEdge[e]:=true IF markedNode[w]==false THEN markedNode[w]:=true S.push(w) ELSE w:=S.pop() a.backtrack(S.top(),w) This algorithm maintains a stack, the top of which is the node currently being explored. The first edge of the current node is taken and followed. This continues until a node is found, which has already been visited (visted nodes are marked). Then the algorithm "backtracks" to the last node in the stack, which still has unmarked edges. Since this procedure only finds nodes reachable from the start node, all nodes are given the chance of being the start node ("FOREACH s in V..."). If a node is found which has not yet been visited by a previous round, it is chosen as a new start node. Note that S.top() has to be nil if the stack is empty, because backtrack will also be called for the start node. The algorithm runs in time O(m+n), because every edge and every node is explored exactly once. If this algorithm is to be run on an undirected graph (represented as a bidirected graph), then every marking of an edge also has to mark the reverse edge. The marking of the visited edges can be avoided if adjacency arrays are used: For each node, one stores not only the index of its edges in the edge array, but also the index of the last edge visited. One visits the edges one after the other and increases the last-visited-edge-pointer each time. Thus, each node knows which edge to visit next without the need to mark the edges. New node of a graph traversal: A node which has not yet been visited. In DFS, these nodes are unmarked. Open node of a graph traversal: A node which has been visited and will be checked out again for successor nodes. In DFS, these nodes are marked and kept on the stack. Finished node, Closed node of a graph traversal: A node which has been visited and will not be checked again. In DFS, these nodes are marked and are no longer on the stack. DFS-tree of a graph (V,E): The graph (V,E'), where E' is the subset of those edges over which DFS backtracks. E' can be obtained by plugging the following DFSAction into the DFS-Algorithm: CLASS DFSTree EXTENDS DFSAction parent : Array[n] of Node PROCEDURE root(s : Node) parent[s] := s PROCEDURE backtrack((v,w) : Edge) parent[w] := v When an edge (u,v) is backtracked, then the node u stores a parent-pointer to the node v. The set E' can now be derived from the parent-pointers: E' := {(u,v) | parent[v]==u} DFS-Numbering: A DFS, which assigns to each node a number indicating the time when it is explored. DFS-Numbering can be done by plugging the following DFSAction into DFS: CLASS DFSNumbering EXTENDS DFSAction dfsNum : Array[n] of Integer := <-1,...> d : Integer PROCEDURE root(s : Node) dfsNum[s] := 0 d := 1 PROCEDURE traverse((u,v) : Edge) IF dfsNum[v]==-1 THEN dfsNum[v] := d d := d + 1 Whenever an edge (u,v) is traversed and has not yet been seen, then the DFS-Number d is assigned to v and d is increased. Tree edge of a DFS-run: An edge leading to a new node. That is: For the edge (u,v), v is unmarked. Tree edges make up a spanning tree. Spanning tree of a connected graph (V,E): A tree (V,E'), where E'=dfsNum[v]. Back edges lead from a node to its (grand-)parent in the spanning forest. Forward edge of a DFS-run: An edge leading to a closed node with a greater dfsNum. That is: For the edge (u,v), v is marked and no longer on the stack and dfsNum[u]dfsNum[v]. cross edges connect two nodes, none of which is the (grand-)parent of the other in the spanning forest. A DFS in an undirected graph cannot encounter any cross edges: Since v is marked, it must have been visited. Since dfsNum[u]>dfsNum[v], u must have been unvisited at this time. Now all edges of v will be explored, in particular also {v,u}, which exists in the undirected graph as well. Thus, {v,u} becomes a tree edge -- and cannot be a cross edge. Spanning forest by network problem: The problem of, given a set V of n nodes, O(n) space in memory and m edges one after the other, calculating a spanning tree of the resulting graph. The problem is that the edges cannot just all be read and inserted into the graph, because one only has O(n) space. The problem can be solved as follows: One reads a set E[1] of 2*n edges. Then, one calculates a spanning forest of (V,E[1]) by DFS. This takes time O(n+m), i.e time O(n + n)=O(n) here. It results in the tree edges F[1], where |F[1]| < |V| < |E[1]|. Now, one reads the next set of edges E[2], and computes again a spanning forest of (V,F[1] \/ E[2]). This is repeated until all edges have been read. If there are m edges, then this means repeating an O(n) process m/n times. The resulting total running time is O(n+m). Topological Sorting of a graph: The construction of a list of all nodes, such that any edge of the graph goes from a prior node to a later node. To sort a graph topologically, it must not contain any cycles (it must be a DAG). A topological sorting can be achieved by plugging the following DFSAction into the DFS: CLASS TopSort EXTENDS DFSAction topologicalOrder : List of Node := nil PROCEDURE traverse((u,v) : Edge) IF (u,v) is a back edge THEN THROW exception("Cycle found") PROCEDURE backtrack((u,v) : Edge) topologicalOrder.pushFront(v) Whenever an edge is backtracked, its end node is pushed on the stack. Thus, only closed nodes are put on the stack. The stack is in topological order, because whenever a node u is pushed on it, there is never an edge from another node v on the stack to u. Proof: Assume there was such an edge * If it was a tree edge, then v would not be closed and would not be on the stack * If it was a back edge, then the algorithm would have signalled a cycle * If it was a forward edge, then v would not be closed and would not be on the stack * If it was a cross edge, then u would already be closed Connected component of a graph G: A connected, vertex induced subgraph of G. That is: A subset of the nodes of G, together with all their edges among each other, which is connected. Intuitively, a connected component is a number of nodes, which are all somehow connected by edges of any directions. SCC, Strongly connected component of a graph G: A maximal, strongly connected, vertex induced subgraph of G. Thus, a SCC is a subset of the nodes of G, together with all their edges. Within this subgraph, there is a path from every node to every other node. The SCC can be determined by plugging the following DFSAction into DFS-Numbering: CLASS SCC EXTENDS DFSNumbering openNodes : Stack of Node := nil openComponents : Stack of Node := nil // These are nodes, each of which represents a component // A component is on this stack, if there may still be new // elements for it component : Array[n] of Node := // This array stores for each node one component representative PROCEDURE root(s : Node) super.root(s) openNodes.push(s) openComponents.push(s) PROCEDURE traverse( (u,v) : Edge) super.traverse((u,v)) IF (u,v) is a tree edge THEN // Every new node is seen as a component per se openNodes.push(v) openComponents.push(v) ELSE IF v in openNodes THEN // v is already open, we have a back edge and a cycle. // Mark all nodes on the cycle as belonging to the // same component WHILE dfsNum[v] < dfsNum[openComponents.top()] DO openComponents.pop() PROCEDURE backtrack( (u,v) : Edge) IF v == openComponents.top() THEN // Backtracking over a component representative // Assign this representative to all open nodes openComponents.pop() REPEAT u:=openNodes.pop() component[u] := v UNTIL u == v This algorithm works due to the following SCC-Cycle-Lemma. SCC-Cycle-Lemma: The fact that two nodes are in the same SCC iff they are on a cycle. Necessity: If two nodes u and v are in the same SCC, then there is a path from every node in the component to every other node in the component. In particular, there is a path from u to v and a path from v to u, which form a cycle together. Sufficiency: If two nodes lie on a cycle, then there is a path from every one of these two to the other one. As a result, they form a strongly connected component. Topologically sorted SCC-algorithm: A modified SCC-algorithm, which uses numbers as component representatives. The numbers are assigned in such a way that if each SCC is collapsed to one node, the numbers indicate a topological sorted sequence of the nodes. The backtrack- procedure of the SCC-algorithm needs to be replaced as follows: CLASS TopSortSCC EXTENDS SCC component_counter := 0 PROCEDURE backtrack( (u,v) : Edge) super.backtrack((u,v)) IF v == openComponents.top() THEN openComponents.pop() REPEAT u:=openNodes.pop() component[u] := component_counter UNTIL u == v component_counter := component_counter + 1 This procedure assigns a component number whenever a component is closed. Since the component numbers are assigned in the wrong order, they need to be reversed after the DFS finished: PROCEDURE renumber() FOREACH v in V DO component[v] := component_counter - component[v] Two-edge-connected nodes in an undirected graph: Two nodes, which lie on a cycle. Two-edge-connected component of an undirected graph: A set of two-edge- connected nodes. This means: If one cuts one edge of a two-edge- connected component, all of its nodes still remain connected. Since SCCs in undirected graphs are two-edge-connected components, the SCC algorithm can be used to find maximal two-edge-connected components. Biconnected nodes in an undirected graph: Two nodes, which lie on a simple cycle. Biconnectedness implies two-edge-connectedness. Biconnected component of an undirected graph: A set of edges of the graph, which lie on a simple cycle. The nodes of any biconnected component are biconnected and also two-edge-connected. That is: If one takes out one of the nodes of a component, the others still remain connected. Biconnected components can be found by plugging the following DFSAction into the DFS-Numbering algorithm: CLASS Biconnected EXTENDS DFSNumbering openComponents : Stack of Edge := nil openEdges : Stack of Edge := nil components : Array[m] of Edge := PROCEDURE traverse((u,v) : Edge) super.traverse((u,v)) openEdges.push({u,v}) IF {u,v} is a tree edge THEN openComponents.push({u,v}) ELSE IF {u,v} is a back edge THEN WHILE dfsNum[v] < dfsNum[openComponents.top()[1]] AND dfsNum[v] < dfsNum[openComponents.top()[2]] DO openComponents.pop() PROCEDURE backtrack((u,v) : Edge) super.backtrack((u,v)) IF {u,v} == openComponent.top() THEN openComponents.pop() REPEAT e := openEdges.pop() component[number(e)] := {u,v} UNTIL e == {u,v} This algorithm works very similar to the SCC-algorithm: Both depend on cycles: Each tree edge creates a new component. Whenever a backedge (a cycle) is found, the involved items are merged to one component. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Shortest paths ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SSSP Problem, Single Source Shortest Path Problem: The problem of, given a start node in a graph, finding all shortest paths from this start node to all other nodes. SSSP-Tree-Lemma: The fact that, if there is a graph and a start node s, such that all nodes are reachable from s and such that the graph does not contain any negative cycles, then there is a spanning tree, all paths of which are shortest paths. Proof: Suppose we mark all edges in the graph, which contribute to a shortest path from s to any other node. This shortest path always exists for each node because each node can be reached from s. So these edges will reach all of the nodes of the graph. Now suppose these edges contain a cycle s ~~> u ~~> v -> u. Then v -> u cannot contribute to any shortest path because u can be reached from s by the much shorter path s~~>u. Hence, the solution to a SSSP problem can be given by a spanning tree. This spanning tree can be given by determining for each node a parent node. Negative cycle in a weighted graph: A cycle, in which the sum of the edge costs is negative. Shortest path length for two nodes u and v: The value * my(u,v) := +oo if there is no path from u to v * my(u,v) := -oo if there is a path from u to v which contains a negative cycle * my(u,v) := min({c(p) : p is a simple path from u to v}) else It makes sense to assign my(u,v)=-oo for negative cycles, because if u~~>w~~>w~~>v is a path with the negative cycle w~~>w, then this cycle can be passed infinitely often, lowering the path costs to -oo. The value my(u,v) will be used to argue about algorithms for the SSSP problem, but it will not be used in the algorithms themselves, because the value of my(u,v) is usually unknown. Minimal distance of a node v: The length of the shortest path from a given start node s to v, i.e. my(s,v). Generic SSSP algorithm: The following blueprint for an algorithm which solves a SSSP problem: FUNCTION genericSSSP((V,E,c) : WeightedGraph, s : Node) : (Array of Real,Array of Node) d : Array[|V|] of Real := <+oo,...> // d[i] = Tentative distance from Node s to Node i, this distance is always greater than the length of the shortest path from s to i d[s] := 0 parent : Array[|V|] of Node := // Stores for each node a pointer to the parent WHILE Ex v in V : d(v) != my(s,v) DO // This condition is a theoretic one, it cannot yet be checked Choose e in E relax(e) RETURN (d, parent) FUNCTION relax((u,v) : Edge) : boolean IF d(u)+c(u,v) < d(v) THEN // We can reach u with cost d(u), we can pass // (u,v) with cost c(u,v), hence v can be reached with // cost d(u)+c(u,v). Store u as a new parent of v parent[v] := u d[v] := d(u)+c(u,v) RETURN true RETURN false This algorithm uses an array d of tentative distances. d[i]=x means that a path of length x from s to node i has been found. The algorithm continuously "relaxes edges" until every node has its optimal distance. Relaxing an edge means: Check whether the edge is useful for reaching the end node of the edge. This is the case if the distance of the start node plus the cost of the edge is smaller than the tentative distance of the end node. Then, the start node is saved as a predecessor of the end node and d[i] is updated. Relaxation-Lemma: The fact that in the genericSSSP procedure, the tentative distance d[i] of a node i is its minimal distance, if the edges of the shortest path from the start node to node i are relaxed in the order given by the path. This lemma is quite intuitive: We pick the shortest path from the start node to node i. Then we follow this path and sum up the costs of the edges. This is what the algorithm does, it always adds the cost of edge (u,v) to the edge sum calculated so far in d[u] and stores the result in d[v]. Then d[i] will contain the sum of all the edges of the shortest path, and hence nothing else than the minimal distance of i. Parent-Cycle-Lemma: The fact that, if a node v lies on a cycle of the parent-pointers of the genericSSSP procedure, then the minimal distance of v is -oo. Proof: Suppose the cycle of parent-pointers is s<~~u<~~v<~~u. Then these parent pointers have been installed because c(s~~>u~~>v~~>u) < c(s~~>u), i.e. u can be reached with less costs if one travels from s to u, from u to v and from v to u. Now it is c(s~~>u~~>v~~>u) < c(s~~>u) <=> c(s~~>u) + c(u~~>v~~>u) < c(s~~>u) <=> c(u~~>v~~>u) < 0 So u~~>v~~>u is a negative cycle and my(s,v)=-oo. Bellman-Ford-Algorithm: The genericSSSP procedure, where the main loop has been replaced by: FOR i := 1 to |V|-1 DO FOREACH e in E DO relax(e) // Check for edges which can still be relaxed relaxable : Set of Node := {} FOREACH (u,v) in E DO IF relax(u,v) THEN d[v]:=-oo relaxable := relaxable \/ {v} // Spread the distance -oo to all nodes reachable from WHILE relaxable != {} DO u := pickAny(relaxable) FOREACH (u,v) in E with d[v]!=-oo DO d[v]:=-oo relaxable := relaxable \/ {v} This algorithm relaxes |V|-1 times all edges of the graph. This suffices to find the shortest path to each node: The shortest path to any node i has at most |V|-1 elements. In the first loop execution, its first edge will be relaxed (among others). In the second loop execution, the first edge will not change any more, but the second edge will be relaxed (among others). By help of the relaxation lemma, the shortest path is found. After having done the relaxation, those edges which could still be relaxed once more deserve attention: They lie on a negative cycle. This algorithm works for arbitrary edge weights. It could be improved by relaxing only those edges, the start node of which received a different d[i]. It runs in time O(|V|*|E|). Some-to-Some-Problem: The problem of, given a weighted graph (V,E,c) and two subsets S < V and T < V, finding a shortest path with the property that it starts in a node in S and ends in a node in T. This problem can be solved by adding an artificial source node s, which is connected by a 0-weigh edge to all nodes in S. Furthermore, an artificial node t is added and all nodes in T are connected to it. Then, Bellman-Ford can be run on ((V,E,c),s). The shortest path from s to t (as found by Bellmann-Ford) is a path with the desired property. Arbitrage-problem: The problem of, given n currencies with their exchange-rates, finding a tour of currencies, such that exchanging an amount of money to the first currency, then the result to the second etc. and finally back to the initial currency will result in an increase of money. This problem can be interpreted as a graph problem: Each currency is a node and each possibility of exchange between two currencies is an edge. The edges are weighted by the exchange rate. The question is: Does there exist a cycle in the graph, such that the product of the edge weights is greater than 1? Since the known graph algorithms are designed to minimize the sum of the edge weights instead of maximizing their product, we use the following fact: x[1] * x[2] * ... * x[k] > 1 <=> 1/x[1]*1/x[2]* ... *1/x[k] < 1 <=> log(1/x[1]*1/x[2]* ... *1/x[k]) < log(1) <=> log(1/x[1]) + ... log(1/x[k]) < 0 If we replace the edge weights by log(1/x[i]), then the problem reduces to finding a negative cycle. Bellman-Ford can be used to identify such negative cycles. Solution to the SSSP problem on DAGs: The genericSSSP procedure, where the main loop has been replaced by: topologicalOrder := topologicalSort(V,E,s) FOREACH v in topologicalOrder DO FOREACH (v,w) DO relax(v,w) If the nodes are in topological order, it follows from the relaxation lemma, that all shortest paths are found. The runtime is O(|V|+|E|). Dijkstra's Algorithm: The SSSP algorithm, where the main loop has been replaced by: Q : PriorityQueue of Node with Priority d Q.insert(s) WHILE Q != nil DO u := Q.deleteMin() FOREACH (u,v) in E DO IF relax(u,v)==true THEN IF v in Q THEN Q.decreaseKey(v) ELSE Q.insert(v) This algorithm is a kind of greedy BFS: It applies the relax-function first to those nodes, which promise to be part of a short path, i.e. those with a small d[i]. Relax candidates wait in a priority queue, in which those with a small d[i] are in front. The smallest node is taken and its outgoing edges are relaxed. In doing so, the end nodes of these relaxed edges also become relax candidates. They are put into the queue according to their temptative distance. If a node becomes a relax candidate, when it is already a relax candidate, its position in the queue is updated according to the new d[i] (Q.decreaseKey). The first node of the queue does already have its minimal distance d[i]==my(s,i). Proof: Suppose there was a shorter path from s to i. Then some nodes on this path will already have been explored by the algorithm. Be j the last node on this path explored by the algorithm. Since this path is shorter than the one found by the algorithm, d[j] v[1] --0--> v[2] --0--> v[3] ... ...where edges (v[i],v[j]) have to be added for all j>i. These edges have weight c(v[i],v[j])=j-i. If Dijkstra starts exploring this graph from v[0], it will explore all edges (v[0],v[i]) and add all nodes v[i] to the queue, each with an temptative distance of d[v[i]]=i. Then it follows the smallest edge to v[1] and is forced to correct all temptative distances d[v[i]] by 1 etc. Limousine-problem: The problem of, given a map of a city with only east-west or north-south streets, finding a cheap route from one street intersection to another, where the cost of a route is defined as the the length of the route plus the number of turns times a constant C. That is: Find a short path from one point to another, avoiding turns. The problem can be interpreted as a graph problem, with the street intersections being the nodes and the street parts the edges. The edge weights are the lengths of the street parts. In order to weight the turns, we split up the graph into two parts: The graph of all east-west-streets and the graph of all north-south-streets. For each intersection of a north-south-street with a east-west-street, we add an edge of weight C from the one graph to the other. Now, shortest paths can be calculated e.g. by Dijkstra and their cost will automatically include the cost for turns. SSSP-in-SCC-problem: The problem of solving the SSSP in a graph (V,E) with SCCs of at most a constant size C. This problem can be solved in time O(|V|+|E|): First, all SCCs need to be calculated (O(|V|+|E|)). Next, the shortest paths from all nodes within the SCC to all other nodes within the SCC are calculated. Due to the constant number of nodes in the SCC, this takes constant time. Then, the SCCs can be "unwringled": We make a copy A' of all nodes A in a SCC, and delete all edges in the SCC. We only keep the edges entering the SCC (connecting them to their node in A) and those leaving the SCC (connecting them from their node in A'). Now, we add edges between A and A', which mirror the length of the shortest path from each node in A to each node in A'. Thereby, the whole graph becomes acyclic. Now, the SSSP-problem can be solved in time O(|V|+|E|). Left-right-minimum of a sequence s of type Real: A value, which is smaller than all values with a lower index. Example: s=<5,3,8,4,2,7> ^ ^ <-- left-right-minima Probability of a value of being a left-right-minimum: The value P(x[i]==min({x[1],x[2],...x[i]})) = 1/i where x is the sequence and i is the position of the value, counting 1-based. Expectation of the number of left-right-minima: The value P(s[1] is LRmin or s[2] is LRmin or ...) = 1/1 + 1/2 + 1/3 + ... 1/s.length = H(s.length) < ln(s.length)+1 where s is the sequence of type Real. Expectation of the number of decreaseKey operations in Dijkstra's Algorithm on a graph (V,E): A value smaller than |V|*ln(|E|/|V|). Proof: Consider a node v with all its entering edges e[i]. We give numbers to the edges in the order in which Dijkstra explores them: e[1] is the edge explored first etc. Let d[i] be the tentative distance of the start node of e[i]. Then a decreaseKey-operation takes place for every left-right-minimum of the sequence The number of left-right minima is H(k). The first node is always a left-right-minimum, but causes it an insert-operation instead of a decreaseKey-operation. hnce we have to count H(k)-1 decreaseKey- operations per node.The expectation for the number k of edges of a node is |E|/|V|. Since we have |V| nodes, the expectation of the number of decreaseKey-operations is |V|*(H(|E|/|V|)-1) < |V|*ln(|E|/|V|) Discrete function: A function, the range of which is finite. Discrete priority queue: A priority queue, the priority function of which is discrete. Usually, the priority function maps to a finite subset of the natural numbers. Monotone priority queue: A priority queue, in which there is never an element inserted before the element, which has been deleted most recently. That is: "Minima never decrease", the first element of the queue will always have a lower priority value than any element which will take up this first place later. The priority queue of Dijkstra's algorithm is monotone, because the priority values d[i] increase with the distance from the start node. "Range C restricted priority queue": A priority queue, where the difference of the smallest and the largest priority value never exceeds the value C. The priority queue of Dijkstra is restricted to a range of the maximum edge weight: Whenever the first node i is taken out of the queue in order to relax its edges, the new nodes resulting from this relaxation will have a priority value of d[i]+c(e), which is smaller than d[i]+max({c(e) : e in E}). Thus, the maximum difference of the smallest and the greatest node is the maximum edge weight. // I don't know if this has a name, but I feel we need it for // the bucket queue. Bucket queue of type t with priority f: A discrete, monotone and range C restricted priority queue of type t and priority f, which is defined by CLASS BucketQueue of t B : array[C+1] of List of t := minIndex : Integer := 0 PROCEDURE insert(x : t) B[f(x) mod (C+1)].pushFront(x) FUNCTION deleteMin() : t minIndexBound := minIndex + C + 1 WHILE B[minIndex mod (C+1)]==nil DO IF minIndex > minIndexBound THEN RETURN nil minIndex := minIndex + 1 RETURN B[minIndex mod (C+1)].popFront() PROCEDURE decreaseKey(x : t) insert(x) The "" can be done efficiently, if each element stores a reference to its list entry in the bucket queue. The bucket queue keeps C+1 "buckets" (lists) of elements. minIndex points to the bucket, which contains the element with the lowest priority value. The bucket right to B[minIndex] contains the elements, the priority value of which is the minimum priority value + 1 etc. If the end of the array is reached, the sequence wraps to the beginning, i.e. an element x can always be found at position f(x) mod (C+1). If an element is inserted, it is added to the list of its corresponding bucket. If the minimum is to be deleted, minIndex seeks to the right until a new non-empty bucket has been found, which contains the actual lowest priority value element. This element is removed and returned. Since the queue is monotone, new elements will always be inserted to the right of minIndex (modulo C+1). Since the queue is range restricted, the greatest possible difference between the lowest priority value and the largest priority value is C. Hence, a new element with a very large priority value will never overrun the invisible limit of the minIndex-bucket, since its distance to minIndex is C and the queue has a length of C+1. If m is the number of elements in the bucket queue, then it needs a space of O(m+C) and a time of O(C) for insert and decreaseKey. If used in Dijkstra, the time for insert and decreaseKey can be reduced further as follows: Be m the weight of the lightest edge m := min({c(e) | e in E). Determine the ratio of the largest edge weight to the smallest edge weight r := max({c(e) | e in E}) / m Now, we use r+1 buckets of "width" m, that is the first bucket takes nodes v with 0 < d[v] < m, the second bucket takes nodes v with m =< d[v] < 2*m etc. Then the insert-procedure looks as follows: PROCEDURE insert(x : Node) B[d[x]/m mod (r+1)].pushFront(x) ...and the other operations are adjusted appropriatly. This algorithm works because the d[i] have at least a distance of m among each other, since they are sums of edge weights and the minimum edge weight is m. msd, most significant distinguishing bit of two integers x and y: The 0-based index of the left-most bit of the binary representation of x and y, for which x and y are different. Example: msd(53,59): 53 = 00110101 59 = 00111011 ^-- msd, third bit from right The msd can be calculated as follows: msd(x,y) := (x==y) ? -1 : round_down(log2(x XOR y)) If the numbers are in a constant range, then log2(x) can be calculated in O(1). The msd can be seen as a logarithmic distance. Radix heap of type t with priority f: A discrete, monotone and range C restricted priority queue of type t and priority f, which is defined by CLASS RadixHeap of t and priority f K : integer := 1 + round_down(log(C)) B : array[-1,...K] of list of t minVal : integer := 0 PROCEDURE insert(x : t) B[min(K,msd(f(x),minVal))].pushFront(x) PROCEDURE decreaseKey(x : t) insert(x) FUNCTION deleteMin() : t i := min({j : j in {-1,...K} and b[j] != nil }) IF i==-1 THEN RETURN B[-1].popFront() minVal := min({ f(x) : x in B[i]}) minElement := eps({x| f(x)==minVal}) B[i].delete(minElement) FOREACH v in B[i] DO decreaseKey(v) RETURN minElement The radix heap works similar to a bucket queue, but the buckets are "wider" for areas far away from the minimum. Thus, two elements x1 and x2 may fall into two separate buckets, if their priority is close to the element with the lowest priority value. But they may fall into the same bucket, if their priority is far away from the minimum priority. A new element is inserted in function of the logarithmic distance of its priority value to the lowest priority value. A priority value is changed (decreaseKey) by deleting the element and re-inserting it. Deleting it may be implemented efficiently if the element stores a reference to its position. Bucket -1 stores elements, the priorities of which are equal to the lowest priority value (see definition of msd). If the minimum is to be deleted from the queue, the first non-empty list is seeked and its minimum is returned. Then, the queue has to be re-arranged: Each element of the minimum bucket is deleted and re-inserted. Thereby, the old minimum value (the one just removed) is taken as a measure for msd. All elements to the right of bucket B[i] are still in the right order. Proof: Be oldMin the old minimal priority value, according to which the elements are spread over the buckets. Be x the element in bucket B[i] to be removed. Be y an element in the bucket B[i+k]. Then we have to prove, that y still remains in B[i+k], even if x is taken as a measure for msd: i+k = i+k log2(y-oldmin) = log2(x-oldmin) + k log2(y-oldmin)-log2(x-oldmin) = k log2((y-oldmin)/(x-oldmin)) = k (y-oldmin)/(x-oldmin) = 2^k y-oldmin = 2^k*(x-oldmin) y = 2^k*(x-oldmin) + oldmin Now the new position of y is log2(y-x) = log2(2^k*(x-oldmin)+oldmin-x) = log2((2^k-1)*(x-oldmin)) = log2(x-oldmin) + log2(2^k-1) ~= i + k High Part of an integer x with a bit-breadth K: The value highPart(x,K) := x div 2^(K/2) This is the value, which can be obtained by cutting away the right half of the binary representation of x. Low Part of an integer x with a bit-breadth K: The value lowPart(x,K) := x mod 2^(K/2) This is the value, which can be obtained by cutting away the left half of the binary representation of x. Van Emde Boas Search Tree, VEB-Tree of bit-breadth K: A type, which allows the storing and retrieval of natural numbers. The log2 of the numbers may not exceed K. K must be a power of 2. The VEB-Tree is either * the value nil OR // VEB-tree is empty * an integer value OR // VEB-tree just stores 1 value CLASS VEB minM : integer // stores minimum value of values contained maxM : integer // stores maximum value highBits : VEB-Tree of bit-breadth K/2 lowBits : Array[2^(K/2)] of VEB-Tree of bit-breadth K/2 FUNCTION minVal() : integer // Accessor function for minimal value IF this==nil THEN RETURN +oo IF this== e THEN RETURN e RETURN minM FUNCTION maxVal() : integer // Accessor function for maximal value IF this==nil THEN RETURN -oo IF this== e THEN RETURN e RETURN maxM PROCEDURE insert(x : integer) IF this==nil THEN // make trivial VEB-tree of 1 element this:=x RETURN IF this== e THEN // convert trivial VEB-tree to complex one minM := e maxM := e highBits.insert(highPart(e,K)) lowBits[highPart(e,K)].insert(lowPart(e,K)) // Insert x IF lowBits[highPart(x,K)] is empty THEN highBits.insert(highPart(x,K)) lowBits[highPart(x,K)].insert(lowPart(x,K)) IF xmaxM THEN maxM := x PROCEDURE remove(x : integer) IF this==x THEN this:=nil RETURN subtree := lowBits[highPart(x,K)] subtree.remove(lowPart(x,K)) IF subtree is empty THEN highBits.remove(highPart(x,K)) // Did we destroy our own maximum or minimum? IF x==minM THEN minM := lowBits[highBits.minVal()].minVal() IF x==maxM THEN maxM := lowBits[highBits.maxVal()].maxVal() // Check out whether we are trivial IF highBits== x1 AND lowBits[x1]== x2 THEN this := x1*2^K+x2 PROCEDURE locate(x : integer) : {0,...2^K-1} // Returns x if x is in the tree // Otherwise, return next greater element which is in the tree // If there is none, return +oo IF this=={} THEN RETURN +oo IF this== x THEN RETURN x IF x>maxVal() THEN RETURN +oo IF x F : Set of Node := {} // This set will keep the nodes, which already have their final distance Q.insert(s) WHILE Q != nil DO WHILE F != {} DO u := pickAny(F) FOREACH (u,v) in E DO IF relax(u,v)==true THEN IF v in Q THEN Q.decreaseKey(v) ELSE Q.insert(v) F := F \/ deleteSettled(Q) FUNCTION deleteSettled(Q : RadixHeap of type t with priority f) : Set of t // Deletes and returns all elements of Q having a minimal // priority value F' : Set of t := {} i := min({j : j in {-1,...K} and Q.B[j] != nil }) IF i==-1 THEN F' := { x : x in Q.B[-1] } Q.B[-1] := nil RETURN F' minVal := min({ f(x) : x in B[i]}) // Run through all nodes in B[i], select the one with d[v]=minVal // as well as all those nodes which also already have their final // distance FOREACH v in Q.B[i] DO IF f(v) < minVal + minIn(v) THEN F' := F' \/ {v} ELSE decreaseKey(v) RETURN F' This variant of Dijkstra always reloads a couple of nodes from the queue, which already have their final distance. Then the algorithm works on these interesting nodes (kept in F) and afterwards refills F. The algorithm knows that the minimal node v in the queue already has its final distance. If there is another node w after v in the queue, then the algorithm looks at the ingoing edges of w: If the smallest of them comes from a node u, which is before v (i.e. d[u]=d[w]-smallestEdge < d[v]), then this node u already had its final distance, hence also w must have its final distance. Thus, w can be taken into F. By this "caching", the algorithm avoids unnecessary queue-reordering-operations and runs in time O(m+n). Proof: We have to count the expected number of moves of one node in the queue. Be C the maximum edge weight, be K the length of the bucket queue, i.e. K:=1+log(C). Any node v will only move down to minIn(v), because then it is caught and put to F. In the radix heap, this means that the buckets the node will pass on its way from bucket K+1 to its exit position is K + 1 - log(minIn(v)). Assume that log(c(e)) is always rounded down: E(#moves) =< SUM {E(K + 1 - log(minIn(v))) | v in V } =< n + SUM {E(K - log(c(e))) | e in Edges} Let's look at one edge e. The term E(K - log(c(e))) is P(log(c(e))==K-1)*1 + P(log(c(e))==K-2)*2 ... + P(log(c(e))==0)*K = P(c(e) in [2^(K-1),2^K-1] )*1 + P(c(e) in [2^(K-2),2^(K-1)-1] )*2 ... + P(c(e) in [0,1] )*(K-1) = 1/2 * 1 + 1/4 * 2 ... + 1/2^K * K =< 2 Thus, the big term reduces to =< n + m*2 // This proof has been invented by Thomas Schultz, thanx! APSP Problem, All Pairs Shortest Paths Problem: The problem of, given a graph, finding the shortest path from each node to each other node. This problem could be solved by solving the SSSP problem for every node. In the case of negative edge-weights, Bellman-Ford is required. If one would like to use the faster Dijkstra, one has to change the negative edges first. This change may not destroy eventual shortest paths. This is done by node-potentials. Node-potential in a weighted graph (V,E,c): A function pot:V->Real, such that for all pairs of nodes u,v pot(u)+c(u,v)-pot(v) >= 0 This corresponds to "lifting up" all edge weights in order to achive non-negative edge weights. Reduced cost in a weighed graph (V,E,c) with a node-potential pot: The function tau:E->Real, such that tau(u,v) = pot(u)+c(u,v)-pot(v) This is a new cost function, which follows the "lifting up" of the node-potential. Reduced cost of a path p in a weightd graph with node-potntial: The sum of the reduced costs of all edges of p, noted tau(p). Reduced Cost Lemma: The fact that shortest paths in a weighted graph (V,E,c) are the same as in the weighted graph (V,E,tau), where tau is the corresponding reduced cost function. Proof: tau(p) = tau(p[1],p[2]) + tau(p[2],p[3]) + ... = pot(p[1])+c(p[1],p[2])-pot([p2]) + pot(p[2])+c(p[2],p[3])-pot([p3]) + ... = pot(p[1])+ c(p) -pot(p[p.length-1]) Since for a shortest path between two nodes, pot(p[1]) and pot(p[p.length-1]) are constant, shortest paths with respect to tau are exactly the shortest paths with respect to c. Reduced Cost Node Potential Lemma: The following fact: If a weighted graph (V,E,c) does not contain any negative cycles and if s is a node, from which all other nodes can be reached and the node potential is chosen pot(v) := my(s,v), then tau(e)>0 for all edges e. Proof: For any edge (u,v), it holds that my(s,u)+c(u,v) >= my(s,v) <=> my(s,u)+c(u,v)-my(s,v) >= 0 <=> pot(u) +c(u,v)-pot(v) >= 0 <=> tau(u,v) >= 0 APSP finder: The following algorithm, which solves the APSP problem: FUNCTION APSPfinder((V,E,c) : WeightedGraph) : Array of Array of Node // Returns a 2-dimensional array, in which at position i,j // there is the node to go to if one wants to go from // node i to node j Result : Array[|V|][|V|] of Node // Introduce a new node s V := V \/ {s} // Introduce edges from s to all nodes E := E \/ {(s,v) | v in V\{s}} // These edges have weight 0 c := c \/ {((s,v),0) | (s,v) in E} // Run Belly-Ford to get minimal distances from s // to all nodes d := BellmanFord((V,E,c),s).d // Introduce the node potential to be the minimal // distance to s pot := {(v,d[v]) | v in V} // Define tau as usual tau := {((u,v),pot(u)+c(u,v)-pot(v)) | (u,v) in E} // Run Dijkstra for each node FOREACH v in V DO Result[v] := Dijkstra((V,E,tau),v).parent RETURN Result This procedure first finds some node potential function. Then it defines a reduced cost function based on this node potential. Now, Dijkstra can run, because no negative edges are left. Dijkstra runs for every node, thus producing a total running time of O(n*m + n*(m+n*log(n))) = O(n^2*log(n)). Heuristic in a weighted graph (V,E,c) for a target node t: A function f:V->Real, such that f(v)==0, but this is not mentioned in the slides // (?) A*: The following algorithm, which finds a shortest path from one node to another in a weighted graph: FUNCTION A*((V,E,c) : WeightedGraph, s : Node, t: Node, f : Heuristic for (V,E,c) and t) : Sequence of Node d : Array[|V|] of Real := <+oo,...> d[s] := 0 parent : Array[|V|] of Node := Q : PriorityQueue of Node with Priority d(v)+f(v) := <> Q.insert(s) WHILE !Q.empty() DO u := Q.deleteMin() FOREACH (u,v) in E DO IF relax(u,v)==true THEN IF v in Q THEN Q.decreaseKey(v) ELSE Q.insert(v) RETURN FUNCTION relax((u,v) : Edge) : boolean IF d(u)+c(u,v) < d(v) THEN parent[v] := u d[v] := d(u)+c(u,v) RETURN true RETURN false This algorithm picks those nodes first, which not only have a close distance to s, but also promise to have a short distance to t. This algorithm is equivalent to Dijkstra with a node potential of pot(v) = -f(v). It will only work if tau(e)>=0. This is the case if tau(u,v) >= 0 <=> pot(u)+c(u,v)-pot(v) >= 0 <=> -f(u,t)+c(u,v)+f(v) >= 0 <=> c(u,v)+f(v) >= f(u) This is a kind of triangle-inequality. // For another A* approach, see "Intelligence artificielle" // (ia.txt, French) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Maximum Flows ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Capacity graph: The type CLASS capacityGraph EXTENDS Graph capacity : Node -> Real+ sourceNode : Node targetNode : Node This is a weighted connected graph, the cost function of which is called "capacity". The capacity must be positive. Such a graph can be thought of as a network of tubes, where the diameter of the tubes is given by the capacity of the edges. Common Notation: (V,E,cap,s,t) Flow in a capacity graph (V,E,cap,s,t): A function f:E->Real, such that * f(e) >= 0 for all e in E * f(e) =< cap(e) for all e in E * SUM {f(u,v) | (u,v) in E } = SUM {f(v,w) | (v,w) in E} for all v in V\{s,t} The flow can be thought of as water which flows through the tubes. cap(e) gives the maximum number of liters per second, whereas f(e) gives the actual number of liters per second. The last condition states that the amount of water flowing into a node must also leave it again. This does not apply to s and t. Flow graph: A capacity graph with a flow. A flow graph is usually drawn like a normal graph with the notion "cap(e)/f(e)" with each edge. Excess of a node v in a flow graph (V,E,cap,s,t,f): The difference of the flow of its incoming egdes and the flow of its outgoing edges. excess(v) = SUM {f(u,v) | (u,v) in E } - SUM {f(v,w) | (v,w) in E} By definition, the excess of all nodes except for s and t is 0. The excess can be thought of the amount of water running out of a node. ("Ueberlauf"). Value of a flow in a capacity graph (V,E,cap,s,t): The excess of t. Notation: val(f) = excess(t) Maximal flow in a capacity graph: A flow, such that there is no flow with a greater value. The goal of this section is finding a maximal flow. s-t-Excess-Lemma: The fact that the excess of the source node in a capacity graph is the negative excess of the target node. excess(s) = -excess(t) in (V,E,cap,s,t,f) Proof: excess(s) + excess(t) = SUM {excess(v) | v in V } because excess(v)=0 for every v with v!=s, v!=t. SUM {excess(v) | v in V } = 0 because every f(e) appears twice in the SUM: Once as an input to a node and then (with inverse sign) as its output. The only values which do not go away are excess(s) and excess(t) => excess(s) + excess(t) = 0 => excess(s) = -excess(t) Saturated edge in a flow graph (V,E,cap,s,t,f): An edge the capacity of which is its flow. f(e) = cap(e) Cut in a graph (V,E): A subset of the nodes V. Edge leaving a cut S: An edge, the start node of which is in S and the end node of which is outside S. Edge entering a cut S: An edge, the end node of which is in S and the start node of which is outside S. Edge of a cut S: The set of its entering and leaving edges. A cut can also be given by its edges, if the graph is connected. s-t-cut in a capacity graph (V,E,cap,s,t): A cut, which contains s and does not contain t. Capacity of a s-t-cut S: The sum of the capacities of the edges leaving S. Notation: cap(S) = SUM {cap(u,v) | u in S, v in V\S} with the capacity graph being (V,E,cap,s,t) Saturated cut in a flow graph (V,E,cap,s,t,f): An s-t-cut, such that * all leaving edges are saturated, i.e. f(e)=cap(e) * all entering edges have flow 0, i.e. f(e)=0 Cut-capacity-lemma: The fact that the capacity of a cut is always greater or equal to the value of a flow. Proof: val(f) = excess(t) = -excess(s) = -SUM {excess(u) | u in S } = SUM {f(e) | e leaves S} - SUM {f(e) | e enters S} = SUM {f(e) | e leaves S} - 0 =< SUM {cap(e) | e leaves S} = cap(S) Saturated-cut-lemma: The fact that if there exists a saturated cut in a flow graph, then the value of the flow is the capacity of the cut. Proof: See proof of Cut-capacity-lemma: For a saturated cut, the "=<" is a "=". Saturated-cut-maximal-flow-lemma: The fact that if there exists a saturated cut in a flow graph, then the flow is maximal. Proof: Assume there was a flow with a greater value. Then this value would contradict the cut-capacity lemma. If a flow is maximal, then there exists a saturated cut, but there may also be un-saturated cuts. Residual network of a flow graph (V,E,cap,s,t,f): The capacity graph (V,E',cap',s,t), where E' is {e | e is an unsaturated edge in E } \/ {(v,u) | (u,v) is an edge with flow in E} Thus, E' takes up all unsaturated edges of E and in addition takes up the reverse edges of all flow carrying edges. For edges of the first kind, one sets: cap'(e) = cap(e) - f(e) For edges of the second kind, one sets cap'(v,u) = f(u,v) The residual network displays where there can be additional flow in the flow graph. This can happen in two kinds of edges: If an edge is not saturated, flow can be added. If an edge already carries flow, then this flow can be undone. This is mirrored by the edges of the second kind. Algorithm: FUNCTION residualNetwork((V,E,cap,s,t,f) : FlowGraph) : CapacityGraph result : CapacityGraph := (V,E'={},cap'={},s,t) FOREACH (u,v) in E DO IF cap(u,v)>f(u,v) THEN E' := E' \/ {(u,v)} cap' := cap' \/ {((u,v),cap(u,v)-f(u,v))} IF f(u,v)>0 THEN E' := E' \/ {(v,u)} cap' := cap' \/ {((v,u),f(u,v))} RETURN result Minimal cut in a capacity graph (V,E,cap,s,t): A s-t-cut with a minimal capacity. That is: A subset of V, which * contains s * does not contain t * has a minimal sum of the capacities of the leaving edges Residual-path-lemma: The fact that if there is a path from s to t in a residual network of the flow graph (V,E,cap',s,t,f), then the flow cannot be maximal. Proof: The residual network reflects flow that can be added. If one adds flow along the path, then the overall flow increases. Vice versa, if there is no such path, then no flow can be added and the flow is maximal. Max-Flow-Min-Cut-Theorem: The fact that if there is no path from s to t in the residual network of a flow graph (V,E,cap,s,t,f), then the set of all nodes reachable from s in the residual network constitutes a minimal cut. Proof: Be S the s-t-cut, which contains all nodes which are reachable from s in the residual network Gf. We know that all edges e leaving S in the original graph must have maximal flow f(e)=cap(e), because these edges do not appear in Gf. Furthermore, we know that all edges entering S in the original graph must have flow 0, because their reverse edge does not appear in Gf. As a result, S is saturated. It follows by the saturated-cut- maximal-flow-lemma that val(f)=cap(S). Hence, there cannot be a s-t-cut with a smaller capacity that cap(S), because val(f) has to go through it. It follows that cap(S) is minimal, S is a minimal cut. Capacity of a path in a capacity graph: The minimum of the capacities of the edges in the path. Augmenting path in a flow graph G: A path from the source node to the target node in the residual network of G. This path shows that the total flow can still be increased by increasing the flow along this path. Max-Res-Cap-Path-Lemma, Maximal Augmenting Path Lemma: The fact that an augmenting path with a maximal capacity can be found in time O(m*log(m)) in a residual network, where m is the number of edges. Proof by algorithm: FUNCTION findMaxAugPath((V,E,cap,s,t) : CapacityGraph) : List of Node m := |E| // Sort edges by their capacity, big capacity first := sortBy(E,cap,DECREASING) // Do a binary search to find the minimum set of edges, // which still allows a path from s to t lasti:= m i := m/2 REPEAT // Find a path with edges e[1],...e[i] E' := {e[0], ...e[i]} p := findPath((V,E'),s,t) // If there exists a path, try smaller i // else try greater i tmp:=i IF p==nil THEN i:=i+abs(lasti-i)/2 ELSE i:=i-abs(lasti-i)/2 lasti:=tmp UNTIL abs(lasti-i)=<1 && p!=nil || i>=m RETURN p This algorithm tries to build up a path of the most powerful edges. To do so, it sorts these edges by their capacity and then does a binary search to find out how many of these edges are needed for a maximal augmenting path. The time is in O(n+m)*O(log(m)) = O(m*log(m)+n*log(m)) = O(m*log(m)), since m>=n in a connected capacity graph. Augmentations-Lemma: The fact that in an integer capacity graph with m edges and a maximum flow value of maxVal, O(m+m*log(maxVal/m)) augmentations of an initial 0-flow suffice to find the maximal flow. Proof: (?) Ford-Fulkerson Algorithm: The following algorithm, which finds a maximal flow in a capacity graph: FUNCTION fordFulkerson((V,E,cap,s,t) : CapacityGraph) : FUNCTION E->Real n := |V| m := |E| f : FUNCTION E->Real := {(e,0) | e in E} REPEAT (V,E',cap',s,t) := residualNetwork(V,E,cap,s,t,f) p := findPath((V,E),s,t) IF p==nil THEN // If t is not reachable from s in Gf, we have a // maximal flow RETURN f // Find out how much we can increase incVal := min({cap'(e) | e in p}) // Increase flow along p FOREACH e in p DO f(e) := f(e)+incVal UNTIL 42==0 This algorithm bases on the residual network: If there is a path from s to t in it, then there is a possibility to push flow from s to t. The algorithm finds such paths and increases the flow along them. With integer capacities, the algorithm has a running time of O((m+m*log(maxVal/m))*m*log(m)), where maxVal is the value of the maximum flow. Reason: According to the Maximal Augmenting Path Lemma, a maximal augmenting path can be determined in time O(m*log(m)). This procedure is necessary for every flow augmentation. According to the Augmentations-Lemma, O(m+m*log(maxVal/m)) of these augmentations suffice. Thus, the claimed running time results. With non-integer capacities, the algorithm may run forever. Example for this (?) BFS-depth of a node in a graph (V,E) relative to a start node s: A natural number, which identifies the time when the node was found by a BFS starting from s. Thus, s gets the number 0 and each "sphere" discovered by BFS gets a new number. Algorithm: FUNCTION BFSDepth((V,E) : Graph, s : Node) : Array of Integer // Essentially do a BFS q : Queue of Node := marked : Array[|V|] of Boolean := marked[s]:=true depth : Array[|V|] of Integer := <0,...> WHILE q != nil DO u := q.popFront() marked[u] := true FOREACH (u,v) in E DO IF marked[v]==false THEN q.pushBack(v) marked[v]:=true // The depth of a node is the depth of its predecessor + 1 depth[v]:=depth[u]+1 RETURN depth Reverse BFS Depth: A BFS-Depth-Algorithm carried out on a graph, the edges of which have been reversed. Algorithm: FUNCTION reverseBFSDepth((V,E) : Graph, s : Node) : Array of Integer RETURN BFSDepth((V,{(v,u) | (u,v) in E),s) Layered subgraph of a graph (V,E) relative to a start node s: A graph (V,E'), where E' contains only those edges, which connect a node of BFS-depth i with a node of BFS-depth i+1. Algorithm: FUNCTION layeredSubgraph((V,E) : Graph, s : Node) : Graph depth := BFSDepth((V,E),s) E' := {} FOREACH (u,v) in E DO IF depth[u]+1==depth[v] THEN E' := E' \/ {(u,v)} RETURN (V,E') This algorithm divides the graph into "spheres", which indicate the distance of the nodes to the source node. Any edge, which goes within one sphere or which goes backwards is cancelled. This algorithm does not destroy shortest paths from s to any target node. Blocking Flow in a capacity graph: A flow, which uses the capacities in such a way, that no other flow may go from the source node to the target node. That is: For every path from the source node to the target node, there is a saturated edge in it. Algorithm: FUNCTION blockingFlow((V,E,cap,s,t) : CapacityGraph) : FUNCTION E -> Real p : List of Nodes := f : FUNCTION E->Real := {(e,0) | e in E} REPEAT IF p.last()==t THEN // We found a path from s to t, breakthrough // Calculate maximal increase incVal := min({cap(e)-f(e) | e in p}) // Increase flow along path FOREACH e in p DO f(e) := f(e)+incVal IF f(e)==cap(e) THEN E := E\{e} // Reset path p := ELSE // Extend path if possible IF Ex (p.last(),w) in E THEN p.pushBack(w) ELSE // If there is no more edge, we're done IF p.last()==s THEN RETURN f ELSE // Retreat (Backtrack), delete the edge from // the path and the set of edges E := E\{ (p.get(p.length()-2),p.last()) } p.delete(p.last()) UNTIL 42==0 This algorithm works similar to a DFS. It tries to find a path from s to t. Then it sends a flow through it and deletes the saturated edges. Next, it finds another possible path. In doing so, it simply destroys useless or saturated edges. As a result, no path from s to t with an unsaturated edge is left. The running time results from the number of operations: * Number of breakthroughs =< m, because each breakthrough saturates one edge * Number of retreats =< m, since one edge is removed * Number of extend operations =< Number of retreats + n*Number of breakthroughs, because each extend * is later undone by a retreat OR * is part of a breakthrough-path For the first case, only as many extends as retreats can appear. In the second case, a breakthrough-path may involve at most n extends. Thus, there may be all in all n times as many extends as breakthroughs. Since each extend falls either in the first case or the second, the claimed upper bound results. Thus, the running time is in O(m+m+m+n*m)=O(n*m). If the capacities of the edges are only 1 or 0, then the algorithm even runs in O(m+n), since every breakthrough saturates all edges on the path. Note: A blocking flow is not necessarily a maximal flow. Dinic's Algorithm, Dinitz' Algorithm: The following algorithm, which finds a maximal flow in a flow graph FUNCTION dinic((V,E,cap,s,t) : CapacityGraph) : FUNCTION E->Real f : FUNCTION E->Real := {(e,0) | e in E} REPEAT Gf := residualNetwork(V,E,cap,s,t,f) Lf := layeredSubgraph(Gf,s) // In the slides, a reverse BFS is done! (why(?)) fb := blockingFlow(Lf) f := f + fb UNTIL fb is 0-function RETURN f This algorithm sees a maximal flow as a sum of blocking flows. The blocking flows are calculated in the residual network, because this graph reflects the edges which still allow flow. In calculating the blocking flows, the algorithm does not care for the edges within one "sphere" in the sense of the layered subgraph. These edges do not contribute efficiently to a flow. As soon as the direct edges are saturated, the residual network changes and the formerly discarded edges come back into play. The loop is executed at most n times (proof to follow(?)), and each loop requires the calculation of the residual network (O(m)) and a blocking flow (O(n*m)), resulting in a total time of O(m*n^2). If the capacities are either 0 or 1, then the algorithm runs in O((m+n)*sqrt(m)) (why(?)). Bipartite graph: A graph (V,E), in which there is a subset V' of V, such that there is no edge within V' and no edge within V\V'. This graph reflects the idea of two different kinds of nodes: V' and V\V'. Edges always connect nodes of one kind to nodes of the other kind. Notation: (V' \/ V\V',E) Integer-Flow-Theorem: The fact that in a flow graph with integer capacities, there is a maximal flow with integer flow values. Proof: Use Ford-Fulkerson to construct the maximal flow and increase the flow by 1 in each iteration. Matching in an undirected graph (V,E): A subset M of E, such that (V,M) has maximum degree 1. This corresponds to a pairwise mapping of nodes to other nodes. Observations concerning flows: * A maximal flow can contain cycles. Consider e.g. ,-------1000/1000----------> a t <------1/1------- s <------1000/1000------------' * There always exists a cycle-free maximum flow, because the flow in each cycle can be reduced by the minimum flow of its edges without disturbing the maximal flow. * If all edges in a capacity graph have unique capacities, this does not mean that the maximal flow is unique. Consider e.g. s -----1/1------> a -----2/1----> b -----4/1----> t '-----3/0--------------------> * If one replaces each edge by an undirected edge in a capacity graph, then the maximal flow may change. Consider e.g. s <----1/0------ t * If all edge capacities in a capacity graph are multiplied by a positive constant c, then the minimal cut remains unchanged. Proof: Be cap'(e)=c*cap(e) the new capacity function. Consider the minimal cut S and any other cut S'. It holds cap(S) < cap(S') c * cap(S) < cap(S') * c c * SUM {cap(e) | e leaves S} < SUM {cap(e) | e leaves S'} * c SUM {c*cap(e) | e leaves S} < SUM {c*cap(e) | e leaves S'} SUM {cap'(e) | e leaves S} < SUM {cap'(e) | e leaves S'} cap'(S) < cap'(S') Thus, even in the changed graph, S has a smaller capacity than S' and hence S is still the minimal cut. * If all edge capacities in a capacity graph are increased by a constant c, then the minimal cut may change. Consider e.g. a minimal cut S with leaving edges e1 and e2, each with capacity cap(e1)=cap(e2)=1. Now consider a non-minimal cut S' with the single leaving edge e3 with cap(e3)=3. If we now add c=1000 to each capacity, the new capacity of S is cap(S) = cap(e1)+1000 + cap(e2)+1000 = 2002 whereas the capacity of the formerly non-minimal cut S' becomes cap(S') = cap(e3)+1000 = 1003 * If all edge capacities in a capacity graph are increased by a constant c, then the value of the maximal flow does not necessarily change by a multiple of c. Consider e.g. s -----3/2------> a -----1/1----> b ------1/1-----> t '-----------------------1/1-----> The flow is 2 here. If the capacities are increased by c=2, then the two edges entering t triple their capacity to 6 all in all, but the edge leaving s will only allow a flow of 5. Thus the flow will increase from 2 to 5, which is not a multiple of c. Spamming problem: See your mailbox. In terms of graph theory: The problem of, given a network with one spamming node q and a set S of peaceful users, determining the connections where a spam filter needs to be installed, such that all users are protected and a minimal amount of money is spent. Installing a spam filter between the network nodes u and v costs c(u,v). The communication network can be seen as a capacity-graph (V,E,cap) with q in V, and S < V. We create an artificial target node t, which we connect to all nodes of S with a capacity of +oo. The goal is to find a minimal q-t-cut: If we block all edges leaving this cut, we get a total isolation of q with a minimal effort. The minimal cut can be found by finding a maximum flow, creating the residual network and determining the nodes reachable from q. This cut does not contain any node of S: The edges from S to t will certainly also appear in the residual network, due to their high capacity. Thus, if any node in S was reachable from q, then t would also be reachable from q and the calculated flow cannot be the maximal flow. Most vital arc, MVA of a capacity graph: An edge, the deletion of which causes the largest decrease in the value of the maximal flow. Some observations: * A MVA does not need be an edge with a maximal capacity.Consider e.g. a <-----100000/0----- s -----1/1------> t * A MVA does not need be the one with a maximal flow. Consider e.g. ,-------1000/1000----------> a t <------1/1------- s <------1000/1000------------' * A MVA does not need to be leaving a minimal cut. Consider e.g. s -----1/1------------------> b ------10000/2--------> t '-----1/1----> a --1/1-----> * A graph may contain several MVAs. Consider e.g. s -------1/1----> a -------1/1-------> t Bipartite Matching Problem: The problem of, given a bipartite graph (A \/ B,E), finding a matching. This corresponds to the idea of creating couples (a,b) of nodes of A and nodes of B, such that a and b are connected. The bipartite matching problem can be reduced to the maximum flow problem: One adds a source node s and connects it to all nodes of A. Then one adds a target node t and connects all nodes of B to it. All edges are given a capacity of 1. Then the problem is to find a maximal flow, in which all edge flows are 1 or 0: The 1-edges are to keep and the 0-edges are to throw away. According to the Integer-Flow-Theorem, such a flow does exist. Dinner problem: The problem of, given n families with a[i] members and m tables with b[j] chairs, finding a seating order such that no two families sit together at one table. This problem can be solved similarly to the Bipartite Matching Problem: One considers each family a node and each table a node. One adds a source node s and edges with capacity a[i] to every family i. One adds a target node and egdes with capacity b[i] from every table to t. Then one adds edges with capacity 1 from every family to every table. Now, finding a maximal flow with integer capacities in this graph will highlight some edges between the families and the tables, indicating where the family members have to go. Preflow in a capacity graph (V,E,cap,s,t): A function f:E->R, such that * f(e) >= 0 for all e in E * f(e) =< cap(e) for all e in E * SUM {f(u,v) | (u,v) in E } >= SUM {f(v,w) | (v,w) in E} for all v in V\{s} The constraint of no excess is relaxed: A node may have positive excess. Thus, a preflow can be seen as a trial of a flow. As a preflow is very similar to a flow, a preflow is often treated like a flow, especially concerning the existence of the residual network and the notions of saturation and excess. Preflow graph: A capacity graph with a preflow. Such a graph can be imagined as a system of tubes with balloons in between. If there is water preflow in a tube, the balloon at the end of this tube inflates. It is under pressure and needs to be released. Push-operation in a preflow graph: The operation of increasing the preflow across a given edge. As a precondition, the edge capacity must allow this increase and the start node of the edge must have an excess greater or equal to the value of the added preflow. The following algorithm does this and adjusts the residual network: FUNCTION push((u,v) : Edge, d : Real, (V,E,cap,s,t, var f) : PreflowGraph, (V,var E',s,t,var cap') : ResidualNetwork) ASSERT d>0 ASSERT cap(u,v)-f(u,v)>=d ASSERT excess(u)>d f(u,v) := f(u,v) + d // update excess, if it is stored as an array // excess(u):=excess(u)-d // update residual network (V,E',cap',s,t) cap'(u,v):=cap'(u,v)-d // Eventually delete (u,v) IF cap'(u,v)==0 THEN E':=E'\{(u,v)} // Eventually add (v,u) IF not (v,u) in E' THEN E':=E'\/{(v,u)} // Set value of residual capacity cap'(v,u) := f(u,v) Saturating push-operation: A push-operation, which saturates an edge. Level-function in a residual network (V,E,cap,s,t): A function d:V->Natural, such that * d(s)=|V| * d(t)=0 * in the residual network of the flow graph, there is no edge (u,v) such that d(u)=d(v)+1. Steep edges are not permitted by the definition of d. Downward edge in a residual network with a level function d: An edge (u,v), where d(u)=d(v)+1. Horizontal edge in a residual network with a level function d: An edge (u,v), where d(u)=d(v). Upward egde in a residual network with a level function d: An edge (u,v), where d(u)=d(v)-1. Active node in a preflow graph: A graph with positive excess. Eligible edge in a residual network with a level function d: An edge, the start node of which has a higher level than the end node. Since step edges are forbidden, this corresponds to a downward edge. Relabeling a node v in a residual network with a level function d: Increasing its value d(v). This corresponds to lifting the node in the landscape of the graph. Generic Preflow-Push-Algorithm, PP-Algorithm: The following algorithm, which finds a maximal flow in a capacity graph: FUNCTION genericPreflowPush((V,E,cap,s,t) : CapacityGraph) : FUNCTION E -> Real // Preflow f : FUNCTION E->Real := {(e,0) | e in E} n := |V| m := |E| // Node levels d : FUNCTION V -> N := {(v,0) | v in E} d(t) := 0 d(s) := n // Calculate residual network (V,E',cap',s,t) := residualNetwork(V,E,cap,s,t,f) // Saturate each source edge by a preflow FOREACH (s,u) in E DO push((s,u),cap(s,u) ,(V,E,cap,s,t,f),(V,E',cap',s,t)) // Select active node WHILE Ex v in V\{s,t} : excess(v)>0 DO // Select eligible edge IF Ex (u,v) in E' : d(u)>d(v) THEN // What flow can we push along it? d := min(excess(u),cap'(u,v)) push(e,d ,(V,E,cap,s,t,f),(V,E',cap',s,t)) ELSE // We have excess, but no downward edge --> relabel u d(u) := d(u) + 1 RETURN f This algorithm first puts as much preflow as possible on each edge, which leaves s. Then it selects a node with positive excess and distributes the excess among the outgoing downward edges. If there are no such downward edges, the node is risen, until downward edges appear. Finally, the superfluous flow will flow back to s. The algorithm leaves much space for heuristics: An active node has to be selected in some reasonable manner and eligible edges have to be selected. Breaking the ties in a PP-Algorithm: Selecting the active nodes and the eligible edges. arbitrary preflow-push-algorithm, arbitrary PP-Algorithm: The PP-Algorithm, where the active nodes and the eligible edges are selected randomly. Steep-Edges-Lemma: The fact that the PP-Algorithm does not produce any steep edges. Proof: If there is a node on a lower level than a node v, then v can push its flow into this lower node until the edge is saturated. Then, the edge disappears in the residual network. Active-Node-Path-Lemma: The fact that from every active node in a preflow graph, there is a path in the residual network to the source node. This lemma guarantees that we can always get rid of superfluous preflow by pushing the flow towards the source node. Proof: Be G=(V,E,cap,s,t,f) the preflow graph. Be Gf=(V,E',cap',s,t) the residual network of G. Be v the node to be analyzed. Be S the set of nodes reachable in Gf from v S := {w | Ex path in Gf from v to w} Be T the set of the remaining nodes T := V\S Then the total excess within S can be calculated as the difference of the inflow and the outflow SUM {excess(v) | v in S} = SUM {f(e) | e enters S} - SUM {f(e) | e leaves S} In Gf, there is no outgoing edge of S ~ (Ex (u,v) in E' : u in S && ~(v in S)) Hence, the flow coming into S must be 0: All e in E : e enters S => f(e)=0 Thus, the sum of the excesses in S must be negative SUM {excess(v) | v in S} = 0 - SUM {f(e) | e leaves S} < 0 Since s is the only node, which is allowed to have negative excess, s must be in S. Max-Level-Lemma: The fact that during the PP-Algorithm, no node attains a level above or equal to 2*n, where n is the number of nodes. This lemma works as an upper bound for the overall height of the graph landscape. Proof: Look at a node v. By the Active-Node-Path-Lemma, there is a simple path from v to the source node s in the residual network. Such a path has at most n-1 nodes. By the Steep-edges-Lemma, the difference in height of two neighboring nodes is at most 1. Thus, the height covered by the path is at most n-1. The height of the source node s is always d(s)=n. Since the sum of both is d(s)+(n-1)*1 = n+n-1 = 2*n-1, d(v) cannot overrun 2*n. PP-Partial-Correctness-Lemma: The fact that, if the PP-Algorithm terminates, then its result is a maximal flow. This is a claim of partial correctness. Proof: First of all, the result of the algorithm is a flow, because the algorithm only terminates if the excess of all normal nodes is 0. In order to show that it is maximal, it suffices to show that there is no path from the source node to the target node in the residual network (by the Residual-Path-Lemma). Since d(s)=n and d(t)=0 and since the length of the path is at most n-1, there would have to be steep edges in this path in the residual network. This contradicts the Steep-Edges-Lemma. Max-Relabels-Lemma: The fact that the PP-Algorithm executes at most 2*n^2 relabel operations. Proof: There are n nodes and each of them is relabled at most 2*n times (by the Max-Level-Lemma). Thus, the total number of relabeling operations can be 2*n^2 at most. Search-cost-Lemma: The fact that total time for the search for eligible edges in the PP-algorithm can be reduced to O(m*n). Proof: If each active node maintains a pointer to its edges, the eligible edges can be selected one after the other, avoiding a search. If the node is relabeled, its pointer needs to be reset. No edge e=(u,v) prior to the current one may ever become eligible again: This could only happen if flow was sent back to the node via the reverse edge (v,u). Since flow only goes downward, (v,u) would be downward, (u,v) would be upward and hence not eligible. Since a node is at most relabeled 2*n times, the total time needed for searching eligible edges amounts to SUM {2*n*degree(v) | v in V} = 4*n*m in O(n*m). Max-Sat-Push-Lemma: The fact that the PP-Algorithm executes at most n*m saturated pushes. Proof: A saturating push removes an edge (u,v) from the residual network. In order to execute a second saturating push, the edge has to give away flow first. This can be done by pushing flow across (v,u). In order to push flow across (v,u), v has to be lifted above u. By the Max-Relables-Lemma, this vice versa overrunning can only happen n times. Since there are m edges, there can be at most n*m saturating pushes. Max-Nonsat-Push-Lemma: The fact that the PP-Algorithm executes at most (1+m)*2*n^2 non-saturating pushes. Proof: Use a "potential function" phi := SUM {d(v) | v is active} Initially, phi is 0. A relabel operation increases phi by 1. By the Max-Relables-Lamma, there are less than 2*n^2 relabels. A saturating push increases phi by at most 2*n, because one new node becomes active and its level is at most 2*n. By the Max-Sat-Push-Lemma, there are at most n*m saturating pushes. Hence, the total increase of phi is at most 2*n^2*1 + 2*n*n*m A nonsaturating push decreases phi by at least 1, since one node pushes its flow away and becomes in-active. phi will never be negative, because it is the sum of positive values. Hence, the maximal number of non-saturating pushes is (1+m)*2*n^2. Since the non-saturating pushes make up the greatest time factor, the PP-Algorithm runs in time O(n^2*m). PP-Termination-Lemma: The fact that the PP-Algorithm terminates. Proof: In each loop, the algorithm either does a push or a relabel operation. The relabel-operations are bound by the Max-relabels-Lemma and the push-operations are bound by the Max-Sat-Push-Lemma and the Max-Nonsat-Push-Lemma. Hence, the algorithm terminates. PP-Total-Correctness-Lemma: The fact that the PP-Algorithm is totally correct. This follows from the PP-Partial-Correctness-Lemma and the PP-Termination-Lemma. Highest-Level-Preflow-Push-Algorithm, High-Level-PP-Algorithm: A PP-Algorithm, which selects always the active nodes at the highest level. To do so, it maintains a priority queue with priority d, which hands out nodes with a high level first. Number of nodes dominated by a node v in PP-algorithm: The number of nodes, which have a level equal to or less than d(v). dominatedNodes(v) := |{u | d(u) =< d(v), u in V}| Tuning parameter of a PP-algorithm runtime-proof: A real value named K. This value will be chosen to be K := sqrt(m) in the High-Level-PP-Theorem. Scaled Number of nodes dominated by a node v in PP-algorithm: The number of dominated nodes, divided by K: d'(v) := dominatedNodes(v)/K d* of a PP-algorithm: The highest level. d* := max({d(v) | v is active}) phase of a PP-algorithm: The set of all push-operations, which are carried out between two changes of d*. Expensive phase of a PP-algorithm: A phase, which involves more than K non-saturating pushes. Cheap phase of a PP-algorithm: A non-expensive phase. K serves as a parameter for distinguishing expensive and cheap phases. Max-Phases-Lemma: The fact that a PP-algorithm has at most 4*n^2 phases. Proof: Every change of d* starts a new phase. d* can only be increased by Relabel-operations. According to the Max-Relabels-Lemma, there are at most 2*n^2 relabel-operations. Thus, d* is only increased 2*n^2 times -- and hence cannot be decreased more than 2*n^2 times. As a result, there is a total upper bound of 4*n^2 changes. This also bounds the number of non-saturating pushes in cheap phases, because there are at most 4*n^2 phases and every cheap phase has at most K non-saturating pushes. Domination potential function of a PP-algorithm: The function phi := SUM {d'(v) | v is active} Phi-Bound-Lemma: The fact that the value of the Domination potential function phi of a PP-algorithm is always >= 0 and initially =< n^2/K. Proof: Since phi is a sum of positive values, it is greater than or equal to 0. Since initially, the source node dominates over all nodes (n) and all the others (n-1) dominate over the others (n-2), we get a sum of dominated nodes of n+(n-1)*(n-2) =< n^2. Since the domination potential function is this value, divided by K, the claimed bound results. Phi-Increase-Lemma: The fact that each relabel-operation and each saturating push increases the domination potential function phi of a PP-algorithm by at most n/K. Proof: The relabel operation changes phi by increasing the d'(v) of one node. A saturating push also changes phi by the d'(v) of one node, because the flow receiving node may become active and thus change d'. In general, it holds for any node v that d'(v) := dominatedNodes(v)/K =< n/K Thus, the above operations can change phi only by at most n/K. Phi-Decrease-Lemma: The fact that the domination potential function phi of a PP-algorithm is not increased by a non-saturating push operation. Proof: If a push across the edge (u,v) is non-saturating, then u got rid of its excess and does no longer belong to the set of active nodes. By contrast, v may become active. But since v is below u, its d'(v) will be smaller than d'(u). Hence, the potential will not increase, but rather decrease. // It seems to be the wrong way round in the slides. Expensive-Phase-Lemma: The fact that an expensive phase of a High-Level-PP-algorithm decreases the domination potential by at least the number of non-saturating pushes, which is >K. Proof: In the High-Level-PP-Algorithm, the nodes at level d* are chosen first. Since d* does not change during a phase, such a phase will empty the level d*. By the phi-decrease-Lemma, each non-saturating push decreases the potential by at least 1. Since by definition, an expensive phase compromises at least K non-saturating pushes, phi decreases by at least K. High-Level-PP-Theorem: The fact that the High-Level-PP-Algorithm runs in time O(n^2*sqrt(m)). Proof: By the Phi-Bound-Lemma, phi is initially smaller than n^2/K. By the Max-Relabels-Lemma, there are at most 2*n^2 relabel-operations. By the Max-Sat-Pushes-Lemma, there are at most n*m saturating pushes. By the Phi-Increase-Lemma, a relabel-operation or a saturating push increases phi by at most n/K. Hence, phi is always below (2*n^2+n*m)*n/K + n^2/K = (2*n^3+n^2+m*n^2)/K Thus, phi can only be decreased (2*n^3+n^2+m*n^2)/K times. By the Expensive-Phase-Lemma, each expensive phase decreases phi by the number of non-saturating pushes. Thus, there can be only (2*n^3+n^2+m*n^2)/K non-saturating pushes in all expensive phases. By the Max-Phases-Lemma, there are at most 4*n^2*K non-saturating pushes in all cheap phases together. This results in a total sum of non-saturating pushes of (2*n^3+n^2+m*n^2)/K + 4*n^2*K This value can be minimized by chosing K = sqrt(m), resulting in a total number of non-saturating pushes of O(n^2*sqrt(m)) Since the non-saturating pushes make up the greatest time factor of the algorithm, the whole algorithm is in time O(n^2*sqrt(m)). Heuristic improvement of the PP-Algorithm: A modification of the PP-Algorithm, which aims at a better performance. Aggressive Relabeling, Local Relabeling: A heuristic improvement of the PP-Algorithm, which relabels a node v by d(v) := 1 + min({d(w) | (v,w) in E}) Instead of increasing d(v) one by one until d(v) overtakes another d(w), the value is increased immediately in such a way that it overtakes the smallest neighbor. This works like a sequence of normal relabelings and saves unneccessary steps. Global Relabeling: A heuristic improvement of the PP-Algorithm, which sets d to d := reverseBFSDepth(Gf,t) where Gf is the residual network of the graph and t is the target node. This happens at the beginning and then periodically every O(m) phases. Thus, the nodes are given their "optimal" level, which encourages the flow towards the target node. Back-Beam: A heuristic improvement of the PP-Algorithm, which treats nodes with a level greater than n apart. Nodes of this kind cannot be connected to the target node in the residual network, as this would entail steep edges. They may only send their flow back towards the source node. Gap heuristics: A heuristic improvement of the PP-Algorithm, which bases on the observation that if there is an empty level, then the nodes above this level will send their flow back towards the source. If there is a level with no nodes (a "gap"), then there cannot be any edge in Gf from the upper level the lower level (no steep edges). Hence, the only paths from nodes in the upper level go to the source node (at level n). Now, the nodes which still have excess will struggle to rise their level until they can send back the flow to the source. This process can be shortened by the following rule IF Ex i : ~(Ex v in V: d(v)=i) THEN FOREACH v in V: d(v)>i DO d(v):=n Thus, a number of previsible steps is avoided. Cost-Capacity-Graph: The type CLASS costCapacityGraph EXTENDS Graph cost : Edge -> Real capacity : Node -> Real+ supply: Edges->Real, such that SUM {supply(v) | v in V} = 0 This graph can be seen as a transportation network, in which the nodes supply a good (supply(v)>0) or demand that good (supply(v)<0). Meanwhile, the edges have a maximum amount of that good, which they can carry (capacity). They also have costs for every unit of the good crossing the edge. The goal is to find the least expensive flow of goods from the suppliers to the demanders. Supply node in a cost capacity graph: A node with a positive supply value. Demand node in a cost capacity graph: A node with negative supply value. Feasible flow, flow in a cost capacity graph (V,E,c,cap,supply): A function f:E->Real, such that * f(e) >= 0 for all edges e * f(e) =< cap(e) for all edges e * excess(v) = -supply(v) for all nodes v All notions relative to flows in capacity graphs are taken over, particularly "excess" and "saturation". A flow in a cost capacity graph can be found by finding maximal flows: We add a source node s and connect s to all nodes v with positive supply with a capacity of supply(v). Then we add a target node t and connect all nodes v with a negative supply to t with a capacity of -supply(v). A maximal flow in this capacity graph outlines a feasible flow. Cost of a flow f in cost capacity graph (V,E,c,cap,supply): The sum of the edge costs, weighted by the flow: c(f) := SUM {c(e)*f(e) | e in E} Min-Cost-Flow in cost capacity graph: The flow with the least cost. Residual Cost network of a cost capacity graph (V,E,c,cap,supply) with a flow f: The cost capacity graph (V,E',c',cap',supply) with * E' = {e | e is an unsaturated edge in E } \/ {(v,u) | (u,v) is an edge with flow in E} * cap'(u,v) = cap(u,v)-f(u,v) if (u,v) unsaturated in E cap'(v,u) = f(u,v) if (u,v) has flow in E * c'(u,v) = c(u,v) if (u,v) unsaturated in E c'(v,u) = -c(u,v) if (u,v) has flow in E This corresponds to a residual network of a capacity graph, to which the approproiate edge costs have been added. If a flow can be undone, then the edge costs in the residual cost network are negative to indicate a potential decrease in cost. Algorithm: FUNCTION residualCostNetwork((V,E,c,cap,supply) : CostCapacityGraph, f : FUNCTION E->Real) : CostCapacityGraph result : CostCapacityGraph := (V,E'={},cap'={},c'={},supply) FOREACH (u,v) in E DO IF cap(u,v)>f(u,v) THEN E' := E' \/ {(u,v)} cap' := cap' \/ {((u,v),cap(u,v)-f(u,v))} c' := c' \/ {((u,v),c(u,v))} IF f(u,v)>0 THEN E' := E' \/ {(v,u)} cap' := cap' \/ {((v,u),f(u,v))} c' := c' \/ {((v,u),-c(u,v))} RETURN result Min-Cost-Cycle-Lemma: The fact that the existence of a Min-Cost-Flow in a cost capacity graph is equivalent to the absence of negative cycles in the residual cost network. Proof: (?). Cycle-Canceling-Algorithm: The following algorithm, which finds a Min-Cost-Flow in a cost capacity graph: FUNCTION cycleCancel((V,E,c,cap,supply) : CostCapacityGraph) : FUNCTION E->Real f := feasibleFlow(V,E,c,cap,supply) Gf := residualCostNetwork((V,E,c,cap,supply),f) WHILE Ex cycle C in Gf: c(C)<0 DO minVal := min({cap'(e) | e in C}) FOREACH e in C DO f(e):=f(e)+minVal Gf := residualCostNetwork((V,E,c,cap,supply),f) RETURN f This algorithm bases on the Min-Cost-Cycle-Lemma: Any feasible flow is taken and then all negative cycles in it are removed. In doing so, the flow still remains a valid feasible flow. Transportation problem: The problem of, given a cost capacity graph with infinitely large capacities, finding a Min-Cost-Flow. This results in a focus on the costs, ignoring the edge capacities. Minimum-Cost-Bipartite-Matching-Problem: The transportation problem in a bipartite cost capacity graph (A \/ B,E,cap,c,supply), where * supply(a)=1 for all a in A * supply(b)=-1 for all b in B * cap(e)=+oo for all e in E For this problem, the desired min-cost-flow-values may only be integer, i.e. f(e)=0 or f(e)=1 for all edges. This mirrors the idea of suppliers in A, which want to get their goods to the demanders in B with a minimal transportation cost. Weight of a matching in an undirected cost graph (V,E,c): The sum of the weights of the edges. c(M) := SUM {c(e) | e in M} for a matching M Maximum-Weight-Matching: The problem of, given a weighted graph (V,E,c), finding a matching with a maximum weight. Approximate-Weighted-Matching-Algorithm: The following algorithm, which finds a matching in a weighted graph, the weight of which is at least half of the weight of a maximum-weight-matching: FUNCTION appWeightedMatching((V,E,c) : UndirectedCostGraph) : Subset of E M' : List of Edge := <> WHILE E' != {} DO // Pick any node v:= eps(V) // Build up a heavy path greedily WHILE degree(v)>0 DO // Find the edge leaving v with a maximum weight (v,w) := eps({(v,u) | c(v,u)=max({c(v,u)|(v,u) in E})}) M'.pushBack(v,w) V := V\{v} v := w // Now we have a set of edges, try building a matching // of all those at even positions weight1 := SUM {c(M'.get(i*2)) | i=0..M'.length/2} // Try building a matching with all odd edges weight2 := SUM {c(M'.get(i*2+1)) | i=0..M'.length/2} // Return the matching with the larger weight IF weight1>weight2 THEN RETURN {M'.get(i*2) | i=0..M'.length/2} ELSE RETURN {M'.get(i*2+1) | i=0..M'.length/2} This algorithm greedily collects edges to paths. It always choses the heaviest edge to go on. Having taken a node, it deletes the node in order to avoid using it again. If it cannot go on, it restarts at another node. Then it tries the two possible matchings contained in the list: a --> b ==> c --> d ==> e --> f ... The two matchings are the set of "-->" edges and the set of "==>" edges. The algorithm takes the matching, which is more heavy. Let's look at a maximal weight matching M*. Consider any edge e* in M*. Either e* is in M', then c(e*) is in c(M*) as well as in c(M'). Now assume e* is not in M'. Let v be the node of e*, which is deleted first by the algorithm (either the start node or the end node of e*). If the algorithm arrives at v, the edge e* is still alive. Now the algorithm will pick the heaviest edge leaving v. As, according to the assumption, e* is not picked, another edge e' is picked with c(e')>c(e*). Now M* cannot contain e', because M* already contains e* and each node has degree 1 in M*. Hence, M' contains a heavier edge than M*. If this argument is repeated for all edges, it results that the sum of the weights of the edges contained in M' is more than the weight of M*. Hence, the weight of the "more heavy half" of M' is greater than c(M*)/2. Multicasting-Network: A weighted directed acyclic graph (V,E,cap) with a source node s in V and a target node set T < V. All nodes except for the source node have the same out-degree. (?) Rate of a multicasting-network: The capacity of the smallest of the minimal s-t-cuts, where t runs through all the target nodes. h = min{cap(S) | S is a minimal s-t-cut with t from T} Information flow in a multicasting-network (V,E,cap,S,T): A function f:E->{0,1}, which maps each edge to a bit, such that f(u,v) = f(v,w) for all u,v,w in V\T\{s} If f(u,v)=b, one says that node u sends bit b to node v. Multicasting-Problem: The problem of, given integer numbers r and t, creating a multicasting-network with rate r, such that the source node can send a maximal number of bits in the same time, which are all received by the t target nodes. The problem is that each node only has r leaving edges and that two nodes may only be connected by one edge. Thus, intermediate nodes have to be created to multiply the information. Coding in a multicasting-network: The idea that the source node does not send only the bits it should send, but applies some function to them. The target nodes undo this function and get the real bits again. With coding, h bits can be sent at a time, if h is the rate of the network. ___________ Example: __/ | Here, each target node _____a_______ v1 _____ t1 | receives either bit a and / / | b, or one of them and the s ---a xor b---- v2 ===-- t2 | value a xor b. By one bit \_____b_______ / | and the value xor b, the 'v3 ===-- t3 | other bit can be |_______| calculated. Linear Coding: Coding in a multicasting-network, which involves vectors associated to each edge (?). The essential thing seems to be that the coding can be done in polynomial time. (?) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Strings ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Sorting algorithm: An algorithm, which produces the sorted sequence of a given set. Rank of an element x in a set S: The index of x in the sorted sequence of S. Bucket sort: The following sorting algorithm for a set of natural numbers in a fixed range: FUNCTION bucketSort(S : set of Integer) : array of Integer // Buckets count occurences of numbers bucket : array[max(S)] of Integer :=<0,0,...0> FOREACH i in S DO bucket[i] := bucket[i] + 1 // Collect buckets to form the ordered sequence result : array[|S|] of Integer resultIndex := 0 FOR bucketIndex := 0 TO bucket.length()-1 DO WHILE bucket[bucketIndex] > 0 DO result[resultIndex] := bucketIndex bucket[bucketIndex] := bucket[bucketIndex] - 1 RETURN result This algorithm has one bucket for each number which may appear in S. It the goes through all the numbers of the set and throws them into the corresponding buckets (by increasing a counter). Then, it goes through all the buckets and constructs the sorted sequence. This algorithm takes time O(|S|+max(S)). Char, Character: A symbol. Usually, a character is encoded by a natural number. Alphabet: A finite set of characters. Ordered alphabet: An alphabet with a linear order on it. By the help of equality, such an alphabet also posesses a simple order. Constant alphabet: An ordered alphabet of a fixed and finite size. Integer alphabet: The alphabet {1,...K} for some natural number K. Any constant alphabet can be represented by an integer alphabet. When an algorithm uses an integer alphabet, one describes its runtime also in dependence of K. // The paper on the Skew-algorithm defines an integer alphabet // slightly different Sentinel character of an ordered alphabet A: A special character, which is smaller than all characters of A. Thus, the character cannot be in A. "Smaller than" is to be read in the sense of the linear order of A. In an integer alphabet, the sentinel character is usually 0. // This corresponds to the final zero of C-Strings and simplifies // sorting String over an ordered alphabet A: A sequence of elements of A, plus the sentinel character of A. We will assume that A be an integer alphabet with the sentinel character 0. The sentinel character is not noted. Notation: |S| := S.length() S = S[i] := for i>=|S| S.alphabet := A "abc" := <'a','b','c'> S1 + S2 := + S2 // "concatenation" Substring from i to j of a string S: The sub-sequence of S, starting at position i and extending to position j-1, plus the sentinel character. Notation: S[i:j) := S[i:j).length() = j-i // Like for intervals, the ")" is probably supposed to indicate // that the character at position j does not belong to the substring. Prefix of size i of a string S: The substring from 0 to i of S. "i-Suffix of a string S", Suffix at position i in S: The sub-string from i to S.length() of S. Notation: S[i..] Lexicographic order: The following linear order on the set of Strings: S1 < S2 :<=> S1[0] < S1[0] || (S1[0]==S2[0] && S1[1..] < S2[1..]) ) Pattern: A (short) string. Occurence of a pattern P in a string S: A substring of S, which equals P. Distinguishing prefix of a string S in a set of strings: The shortest prefix of S, which is not a prefix of another string in the set. If for a string S, there is no such prefix, then the distinguishing prefix of S is S itself. String sorting problem: The problem of, given a set of strings, producing a sorted sequence of them. A string sorting algorithm only needs to look at the distinguishing prefixes of the strings. Standard sorting algrithm: A sorting algorithm, which can sort any ordered set. For string sorting, it is possible to perform better than Standard sorting algorithms, because one can exploit the internal structure of strings and the recusrive definition of their order. Standard sorting algorithms do O(n*log(n)) element comparisons, where n is the number of input elements. Depending on the size of the strings, each comparison may again contribute a significant time factor. Multikey-Quicksort: The following string sorting algorithm: FUNCTION multiKeyQuickSort(S : Set of String, p : Integer := 0) : array of String ASSERT All strings of S have the same prefix of length p // Trivial cases IF |S|==0 THEN RETURN <> IF |S|==1 THEN RETURN // Choose Pivot-element pivot := eps(S) S1 := {s in S | s[p] < pivot[p]} S2 := {s in S | s[p] == pivot[p]} S3 := {s in S | s[p] > pivot[p]} RETURN multiKeyQuickSort(S1,p) + multiKeyQuickSort(S2,p+1) + multiKeyQuickSort(S3,p) This algorithm get a set of strings S, which already share the first p characters. It chooses one arbitrary string of S, called the "Pivot element". Then it partitions S into 3 subsets by looking at the p-th character: Depending on whether this character is smaller than, equal to or greater than the p-th character of the pivot element, the string goes to S1, S2 or S3. Due to the nature of the lexicographical order, this means that all strings in S1 are smaller than the pivot element, all strings in S2 are greater than the pivot element and nothing can be said about the strings of S2. Then, each group is sorted recursively, where p advances only for the group S2. There are two cases: Let's look first at one string s, which survives each recursive call in the set S1. Assuming a perfect choice of the Pivot-element, this set is halved each time. Thus, recursion depth cannot exceed log(n), where n is the number of strings. Since about O(n) strings face the same fate as s, and each of them requires the comparisons, runtime is at least O(n*log(n)). In the other case, the strings share long common prefixes. Then, all strings go to S2. The algorithm will only stopped after having looked at all characters of the distinguishing prefixes of the strings. Let D denote the sum of the lengths of the distinguishing prefixes, then the runtime is O(D+n*log(n)). It can be shown that for a random choice of the Pivot-element, the same runtime order results. LSD-first Radix-sort, Least significant digit first Radix sort: The following string sorting algorithm: FUNCTION LSDradixSort(S: Set of String) : array of String ASSERT eps(S).alphabet is an integer alphabet alphsize := |eps(S).alphabet| buckets : array[2][alphsize] of List of String <<<>,<>,..>..> FOREACH s in S DO bucket[0][0].pushBack(s) switch := 1 FOR len := max({|s| | s in S}) DOWNTO 0 DO switch := (switch+1) mod 2 FOR i := 0 TO alphsize DO WHILE bucket[switch][i] != <> DO s := bucket[switch][i].popFront() bucket[(switch+1)mod 2][s[len]].pushBack(s) result : array[|S|] of String j := 0 FOR i := 0 TO alphsize DO WHILE bucket[switch][i] != <> DO result[j] := bucket[switch][i].popFront() j := j + 1 RETURN result This procedure does a bucket sort for each position ("digit") of the strings. It starts with the rightmost (the "least significant"), i.e. the last character of the longest string (the others automatically have the 0 there). It puts the strings into the buckets. Then, it goes through all buckets (in the right order) and puts the strings into new buckets according to their but-last character. This continues until all positions have been sorted. The running time in O(N), if N is the total number of characters. This algorithm does not make use of the fact that the distinguishing prefix (which may be much shorter than the string) suffices for sorting. MSDRadixSort, most significat digit radix sort: The following string sorting algorithm: FUNCTION MSDRadixSort(S : List of String, p : Integer := 0) : List of String ASSERT All strings share the first p characters ASSERT eps(S).alphabet is an integer alphabet alphsize := |S.alphabet| IF S.length(),<>,...> FOREACH s in S DO bucket[s[p]].pushFront(s) FOR i:=0 TO alphsize-1 DO bucket[i] := MSDRadixSort(bucket[i],p+1) RETURN bucket[0] + bucket[1] + ... + bucket[alphsize-1] This algorithm does repeated bucket-sorts like LSDRadixSort, but it starts at the first character: It puts the strings into buckets according to their first character and then sorts the buckets recursively. The parameter p indicates the next position to analyze. If the size of the set to sort is smaller than the size of the alphabet, the recursion breaks and multikey-Quicksort is used. This happens in order to avoid that initializing the buckets to <> takes longer than the actual sorting process. Taking this into consideration, the runtime is dominated by the FOREACH-loop. The total number of pushFront()-operations for all strings and through all recursions is D, if D is the sum of the lengths of the distinguishing prefixes. This is because after having bypassed the distinguishing prefixes, all strings are in one-element-sets. Thus, running time is at least O(D). In addition, a break of the recursion causes a call to multikey-Quicksort, which has a runtime of O(n*log(n)), where n is the number of strings. Multikey-Quicksort is only called on alphsize strings and thus, the overall runtime is O(D+n*log(alphsize)). Suffix Sequence of a string S: The sorted sequence of all suffixes of S. Suffix sorting problem: The problem of, given a string S, creating its suffix sequence. Example: S="banana" Suffixes ~~~> Sorted suffixes Pos Pos Rank 0 banana 5 a 0 1 anana 3 ana 1 2 nana 1 anana 2 3 ana 0 banana 3 4 na 4 na 4 5 a 2 nana 5 Suffix array of a String S: An array A of integers, where A[i] is the starting position of the suffix with rank i in the sorted sequence. Notation: SuffixRankToPos[i]=j means: The starting position of the i-th suffix is j Example: S = "banana" SuffixRankToPos = <5,3,1,0,4,2> ^------ The second suffix in the sorted sequence is the 3-suffix S[3..] = "ana" Inverse suffix array of a string S: An array A of integers, where A[i] is the rank of suffix S[i..]. That is: The array can be read as an "annotation" for the string S, which notes for each suffix its rank in the sorted sequence of all suffixes Notation: SuffixPosToRank[i]=j means: The rank of the i-suffix is j Example: S="banana" SuffixPosToRank = 325140 ^-- Suffix "ana" is has rank 1 in the sorted sequence of all suffixes (cf. ex. above) The suffix array can be translated to the inverse suffix array and vice versa by doing PROCDURE swap(A : array of Integer) : array of Integer B : array[A.length] of Integer FOR i:=0 TO A.length-1 DO B[A[i]] := i RETURN B // I find the notions "suffix array" and "inverse suffix array" hard // to distinguish. Therefore, I chose different names. // SuffixRankToPos means: I put the rank as an index, so I get the // starting position as a result. (and vice versa) Lexicographic names for a set S: A function f:S->A so that S[i] f(S[i]) triples2 : array[|S|/3] of String := // Put these triples together, sort them triples1&2Sorted : array of String := LSDradixSort({s | s in triples1 or triples2}) // Create lexicographic names for triples1&2Sorted lex : FUNCTION triples1 \/ triples2 -> Integer := {} lex(triples1&2Sorted[0]):=0 FOR i:=1 TO triples1&2Sorted.length-1 DO IF triples1&2Sorted[i]!=triples1&2Sorted[i-1] THEN // Different from previous => new lexicographic name lex(triples1&2Sorted[i]) := lex(triples1&2Sorted[i-1])+1 ELSE // Same as previous => same lexicographic name lex(triples1&2Sorted[i]) := lex(triples1&2Sorted[i-1]) // Translate triples1 and triples2 to their lexicographic // counterparts triples1Lex : array[|S|/3] of Integer triples2Lex : array[|S|/3] of Integer FOR i:=0 TO |S|/3 DO triples1Lex[i] := lex(triples1[i]) triples2Lex[i] := lex(triples2[i]) // Concatenate triples1Lex + triples2Lex triples1&2Lex : array of Integer := triples1Lex+triples2Lex // Create their suffix sequence triples1&2LexSuffixRankToPos : array of Integer // Do we have all distinct elements? IF lex(triples1&2Sorted[triples1&2Sorted.length-1] == triples1&2Sorted.length-1 THEN // Their lexical number is their SuffixPosToRank, // so do a swap swap to get the SuffixRankToPos triples1&2LexSuffixRankToPos := swap(triples1&2Lex) ELSE // Use recursion triples1&2LexSuffixRankToPos := skew(triples1&2Lex) // Convert this back to the original suffixes suffixes1&2Sorted : array[|S|*2/3] of (Suf:String,Pos:integer) FOR i := 0 TO |S|*2/3-1 DO pos := triples1&2LexSuffixRankToPos[i] IF pos < triples1Lex.length THEN pos := pos*3+1 ELSE pos := (pos-triples1Lex.length)*3+2 suffixes1&2Sorted[i] := (S[pos..],pos) // Now take care of the suffixes starting at positions i, // where i mod 3 = 0: Get them and sort them. The suffix of // each of these suffixes has already been sorted in the previous // rounds, so we only need to look at the first character and // the position of the following suffix in the sorted sequence // (This is neglected here) suffix0 : Set of (String,integer) := {(S[i..],i)|i mod 3 == 0} suffix0Sorted : array of (Suf:String,Pos:integer) := sort(suffix0) // Now merge suffixes1&2Sorted and suffixes0Sorted // This merging can be done efficiently by taking into account // the sub-suffixes, which are already sorted // (This is ignored here) result : array[|S|] of integer i := 0 index0 := 0 index1&2 := 0 WHILE i < |S| DO IF suffix0Sorted[index0] SuffixRankToPos := skew(S) L := 0 FOR i:=0 TO |S|-1 DO // Look at the suffix starting at i // Does it have any predecessor? IF SuffixPosToRank[i]>1 THEN // Find its predecessor in the sorted sequence pred := SuffixRankToPos[SuffixPosToRank[i]-1] // Count shared characters, starting at L (!) WHILE S[i+L]==S[pred+L] DO L := L + 1 lcp[SuffixPosToRank[i]] := L L := L - 1 IF L < 0 THEN L := 0 RETURN lcp The main loop is executed |S| times. In each execution, L may be increased by the inner WHILE-loop and is then decreased by 1. Assume the first WHILE-loop-run increases L to n. As L cannot bypass n, each following WHILE-loop-run can at least increment L by one. Thus, we still have a total time of |S|+|S|, i.e. O(|S|). If another WHILE- loop-run wants to increase L more than once, say k times, then the first WHILE-loop-run has to give up exactly k incrementations. But all in all, the total number of incrementations will not bypass n. Thus, the algorithm needs time O(|S|). L only needs to be decreased by 1 each time due to the following LCP-Continuity-lemma. LCP-Continuity-lemma: The fact that, given a string S and a position i, the i-suffix has more shared characters with its predecessor in the sorted sequence than the (i-1)-suffix minus 2. lcp[SuffixPosToRank[i]] > lcp[SuffixPosToRank[i-1]] - 2 Example: S="banana", i=3 lcp[SuffixPosToRank[3]] > lcp[SuffixPosToRank[2]] - 2 lcp[4] > lcp[5] - 2 1 > 2 - 2 Picture: SuffixSeq Rank lcp A 0 0 Ana 1 = SuffixPosToRank[3] 1 -, ANAna 2 3 | banana 3 0 | 1 > 2 - 2 na 4 0 | NAna 5 = SuffixPosToRank[2] 2 -' "ana" has more shared characters with its predecessor "a" than the suffix "nana" has with its predecessor "na", minus 2. Proof: If the (i-1)-suffix shares less than 1 character with its predecessor, then the lemma follows immediatly, because the left-hand side of the inequality is positive, while the right-hand side is negative. Now let's suppose the (i-1)-suffix shares L>1 characters with its predecessor, i.e L=lcp[SuffixPosToRank[i-1]]>1. In the example, i-1 is 2 and the 2-suffix "nana" shares 2 characters with its predecessor "na". Be j the position of this predecessor in the string, i.e. j=SuffixRankToPos[SuffixPosToRank[i-1]-1], which is 4 in this case. Since the j-suffix and the (i-1) suffix share L characters, the (j+1)-suffix and the i-suffix must share L-1 characters ("a" and "ana" share 2-1=1 characters). Since the j-suffix preceded the (i-1)-suffix in the sequence and they share L characters, the (j+1)-suffix must also precede the i-suffix. In the example, the suffix "a" precedes the suffix "ana", although there may be other suffixes in between. The goal is to show that the i-suffix shares more than L-2 characters with its immediate predecessor. We already showed that the i-suffix shares L-1 characters with the (j+1)-suffix, which is a (possibly far) predecessor. But since the suffixes are sorted, all suffixes between the (j+1)-suffix and the i-suffix also have to share at least those L-1 characters. In particular, the direct predecessor of the i-suffix must share at least L-1 characters with the i-suffix, which is the claim. // I owe thank to David Steurer, who helped me understanding this // lemma. Labelled tree: The type CLASS LabelledTree EXTENDS Graph ASSERT the graph is a tree root : Node label:Edges->X d:Node->Y This is a directed tree with a function label, which equips each edge with a name or symbol and a function d, which equips each node with a name or symbol. Suffix-tree of a string S: A labelled tree, so that * its edges are labelled with substrings of S * there is one leaf for each suffix of S * the concatenation of the edge labels on a path from the root to a leaf gives the suffix of the leaf * labels of edges leading to siblings have a distinct first character * the label of a node is the length of the concatenations of the labels of the edges on the path from the root to the node. The suffix-tree can be constructed by hand by the following (simple) algorithm: * Do a suffix sort of S, taking into account the sentinel character. * Introduce a root node and add all suffixes of S as children * While there is still a suffix x!=<0> which is a prefix of sibling suffixes of x * chop off the trailing 0 * add 0 as a child of x * make all sibling suffixes with prefix x children of x * chop off x of all these new children * Label the edges with the node they lead to and compute d Example: S="banana" Suffix sort Make tree root 0 banana0 ~~~~~> 6 0 ~~~~~> |-> 0 1 anana0 5 a0 |-> a -> 0 2 nana0 3 ana0 | '-> na -> 0 3 ana0 1 anana0 | '-> na0 4 na0 0 banana0 |-> banana0 5 a0 4 na0 |-> na -> 0 6 0 2 nana0 '-> na0 The suffix tree can also be constructed by the following (complicated) algorithm: FUNCTION suffixTree(S : String) : LabelledTree SuffixRankToPos := skew(S) lcp := lcp(S) T : LabelledTree := (V={root},E={},label={}, d={}) d(root) := 0 v := root FOR i := 0 TO |S|-1 DO // Insert the suffix with rank i into the tree WHILE d(parent(v)) > lcp[i] DO v := parent(v) IF d(v) > lcp[i] THEN // Split the edge v' := newNode(V) d(v') := lcp[i] e1 := newEdge(E,parent(v),v') label(e1) := label(parent(v),v)[0..lcp[i]] e2 := newEdge(E,v',v) label(e2) := label(parent(v),v)[lcp[i]+1..] E := E \{(parent(v),v)} v := v' // Now d(v) == lcp[i] v'' := newNode(V) d(v'') := |S|-SuffixRankToPos[i] e := newEdge(E,v,v'') label(e) := S[SuffixRankToPos[i]+lcp[i]..] RETURN T This algorithm inserts the suffixes according to their sorted sequence into the tree. Therefore, it can build the edges from the root to the leaves. Only if a suffix is not a prolongation of its predecessor, the algorithm needs to go up in the tree. If it finds a suitable node (lcp[i]==d[i]), it simply adds the new suffix. If not, it breaks the edge in two: The first half gets the first part of the label, the second gets the second part. In between, a new node is added. To this new node, the suffix can then be added as usual. This algorithm works similar to DFS. The main loop loops n:=|S| times. Since each edge is created once and backtracked at most once, we get an additional time consumption of O(2*n), which does not influence the overall running time of O(n). Multikey Suffix Quicksort: The Multikey Quicksort algoroithm, applied to the suffixes of a string. For the string S with n=|S| and sigma= |S.alphabet|, the following holds: * It is not necessary to make copies of all suffixes, because the algorithm can as well work on pointers to the suffixes in S. * If sigma>n, then one can design S such that the sum D of the lengths of the distinguishing prefixes is in O(n): One just creates S such that it does not contain equal characters. As a result, the algorithm does not need any recursive call. * For two non-overlapping substrings S[i:i+k) and S[j:j+k), the probability that they are equal is sigma^-k. Proof: For two randomly chosen characters, the probability of being equal is 1/sigma. If the strings are to be equal, their first character has to be equal and their second character has to be equal etc. Since the choices of characters are independent, it is P(S[i+1]=S[j+1] | S[i] = S[j] ) = P(S[i+1]=S[j+1]) Thus, the overall probability amounts to the product of all single probabilities: 1/sigma*1/sigma*...1/sigma = sigma^-k. * Even if two strings overlap, the probability of being equal is still sigma^-k, because the strings are independently chosen. * The probability of two equal substrings of length k in S is bounded by n^2*sigma^-k. The substring starting at position 0 could be equal to (n-k-1) other substrings. The substring at position 1 could be equal to (n-k-2) substrings etc.. All in all, (n-k +1)(n-k)/2 < n^2 pairs could be equal. For each of them, the probability is sigma^-k. Thus, the probability of any two equal substrings is less than n^2*sigma^-k. * The probability that there are any two suffixes with a common prefix of length 3*log(n)/log(sigma) or longer is at most 1/n. The bound n^2*sigma^-k for the probability of equal substrings of length k also applies to substrings of length k or longer, because every equal substrings of length >k also have to match within the first k characters. Thus, we can set k=3*log(n)/log(sigma) and obtain n^2*sigma^(-3*log(n)/log(sigma)) = n^2*sigma^(log(n^-3)/log(sigma)) = n^2*n^-3 = 1/n * The average case for a multikey suffix quicksort run is bounded by O(log(n)*n). Two cases can be distinguished: * There exist at least two suffixes, which share a prefix of length at least 3*log(n)/log(sigma). The probability of this case is bounded by 1/n. These two strings will force the algorithm to keep them in the same bucket until their equal parts have been bypassed. The sum D of the lengths of all equal parts is bounded by n^2. Thus, we get the running time of O(D+n*log(n)) = O(n^2+n*log(n)). * If there do not exist any two substrings with an equal prefix of 3*log(n)/log(sigma), then D is bounded by n*3*log(n)/log(sigma). The running time is thus O(D+n*log(n)) = O(n*3*log(n)/log(sigma) + n*log(n)) = O(n*log(n)). The average case running time is the sum of the special case running times, weighted with their probabilities: O( (n^2+n*log(n)) * 1/n + n*log(n) * (1-1/n) ) = O(n*log(n)) The prefix length 3*log(n)/log(sigma) has been chosen arbitrarily: If one choses a smaller value, the first case becomes less probable, but the second case becomes worse. Vice versa, if one choses a greater value, then the second case becomes better, but the first case becomes more probable. Burrows-Wheeler-Transform of a string S: The string BW(S), which is obtained as follows: Create a matrix of size (|S|+1) x (|S|+1) and write S with its sentinel character in the first line. Then rotate this line one to the left with the first character wrapping to the last position and write this into the second line. Continue until the matrix is filled up. Now sort the lines alphabetically. The right-most column of the matrix, read top-down, is the Burrows-Wheeler- Transformation of S. Since the sentinel character is smaller than all other characters, the sorting of the matrix lines boils down to the sorting of the suffixes of S. Hence, BW(S) can also be calculated as the sequence of the characters, which precede the suffixes of the sorted suffix sequence. The algorithm would be: FUNCTION BW(S : String) : array of S.alphabet SA := skew(S) result : array[|S|+1] of S.alphabet result[0] := S[|S|-1] FOR i:= 0 TO |S|-1 DO result[i+1]:=S[ (SA[i]-1+|S|) mod |S| ] RETURN result Matching: The process of finding the position(s) in a string, where a pattern occurs as a substring. Brute-Force-Matching: The following (simple) algorithm for finding a pattern P in a string S: FUNCTION bruteForceMatching(P,S: String) : Set of Integer result : Set of Integer := {} FOR j := 0 TO |S|-|P| DO IF S[j..j+|P|] == P THEN result := result \/ {j} RETURN result This algorithm just checks out every position of S. If P matches with all of its characters, then the position is added to the result-set. The algorithm runs in time O(|P|*|S|). Suffix-Tree-Matching-Algorithm: The following algorithm, which finds a pattern P in a string S by using the suffix tree of S: FUNCTION suffixTreeMatch(P,S : String,(V,E,root,label,d) : SuffixTree) : Set of Integer v := root i := 0 // Index of first character, that has not yet been matched WHILE i < |P| DO IF Ex (v,w) in E : label(v,w) == P[i..i+|label(v,w)|-1] THEN i := i + |label(v,w)| v := w ELSE RETURN {} RETURN { |S| - d(v') | v' is a (grand-)child of parent(v) and v' is a leaf} This algorithm goes down in the suffix tree with edges matching the pattern. When the pattern has been eaten, all child suffixes of the current node are occurences. It has a runtime of O(|P|+n), where n is the number of occurences of P in S. The convenient thing about this algorithm is that once the suffix tree of S has been calculated, different patterns can be searched very efficiently. FSA-Matching-Algorithm: The following algorithm, which uses a FSA to find a pattern P in a string S: FUNCTION FSAMatching(P,S : String) : Set of Integer result : Set of Integer := {} myFSA : FSA := (Q : {0,...|P|}, // States == positions in P I : S.alphabet, // Input will be S s0 : 0, // Start state = 0 F : {|P|}, // End state = length of P d) // Transition function myFSA.d(q,c) := max({ j | j=0...|P|, P[0..j) == P[q-j:i) + }) // with max({}) := 0 result : Set of Integer := {} state := s0 pos := 0 WHILE pos<|S| DO state := myFSA.d(state,S[pos]) pos := pos + 1 IF state in myFSA.F THEN result := result \/ {pos-|P|} RETURN result This algorithm first constructs a FSA for the pattern. The FSA has a state for each position of P. It is fed with the string S. With each character, it remains in the start state. When the first character of P occurs, the FSA advances one state. If the second character occurs, the FSA advances one state and so on, until the final state has been reached, meaning that the pattern occured. If the FSA encounters a character c which disturbs the pattern, it returns to a previous state q'. This state is chosen such that it allows continuing a new start of the pattern, which may lie in the pattern itself. If the current state is q, then the new state is the length of the longest prefix of P, which is a suffix of the part of the pattern read so far (P[0..q]) plus the character just read (c). Example: P = "bonbons" S = "bonbonbons" ^--- The FSA runs till here (q=6) until it encounters a mismatch (c='b'). Now, the FSA will not re-start matching from P[0], but it will continue with matching the second 'b' of the pattern with the third 'b' of the string: It knows that the last three characters of the matched part "bonbon" may be the starting of another match of "bonbons". Here, the FSA would be 0 --b--> 1 --o--> 2 --n--> 3 --b--> 4 --o--> 5 --n--> 6 --s--> 7 / \ ^----------b------' '-b-' with d(q,c) = 0 for all transitions not mentioned and F={7} This FSA needs the transition function to be stored, occupying a space of O(|S.alphabet|*|P|). To calculate d, a time of O(|S.alphabet|*|P|^3) is needed, if this is done naively. The running time of the loop is O(|S|). Although the FSA can the be used to search in different strings for the same P, this algorithm is rather inefficient. Morris-Pratt-Failure-Function of a pattern P: A function F:{1,..|P|}->{0,..|P|}, which returns for each position i=0..|P| the length of the longest prefix of S, which is a suffix of P[0:i). F(i) := max({j : j=0..i-1, P[0:j) == P[i-j:i) }) Here, it is max({}):=0. Often, F is stored as a sequence <0,F(1),F(2),...>. Example: P="bonbon" F= _000123 ^--- The 2 preceding characters of this position ("bo") match with the pattern start. The Morris-Pratt-Failure-Function can be calculated as follows: FUNCTION MorrisPrattFF(P : String) : array of Integer F : array[|P|+1] of integer posInP := 1 prefixLen := 0 F[1] := 0 // Run over P WHILE posInP < |P| DO // Run over a matching with P[0..prefixLen] WHILE posInP < |P| && P[posInP] == P[prefixLen] DO posInP := posInP + 1 prefixLen := prefixLen + 1 F[posInP] := prefixLen // No matching found, try next char IF prefixLen == 0 THEN posInP := posInP + 1 // Partial match found, try to use this as a new beginning ELSE prefixLen := F[prefixLen] RETURN F This algorithm does a matching of the pattern in the pattern itself. It runs over the whole pattern and with each character, tries to match with the start of P. While this succeeds, F is filled appropriately. As soon as the complete match fails, prefixLen jumps back to the longest prefix, which has been found -- and this is given by F. During this algorithm, prefixLen is always smaller than posInP, which allows the assignment prefixLen := F[prefixLen]. Since posInP only increases, the running time is in O(|P|). Morris-Pratt-Matching-Algorithm: The following algorithm, which finds a pattern P in a string S by the Morris-Pratt-Failure-Function of P: FUNCTION MorrisPrattMatching(S,P : String, F : MorrisPrattFF(P)) : Set of Integer result : Set of Integer := {} posInS := 0 matchedChars := 0 // Run over the whole string S WHILE posInS + (|P|-matchedChars) < |S| DO INVARIANT A: P[0:matchedChars) == S[posInS-matchedChars:posInS) INVARIANT B: occurences before posInS-matchedChars are in // Run over a sequence of matching characters WHILE matchedChars < |P| && P[matchedChars]==S[posInS] DO matchedChars := matchedChars + 1 posInS := posInS + 1 // All chars matched? ---> husch husch ins Koerbchen IF matchedChars == |P| THEN result := result \/ {posInS - |P|} // No chars matched at all? ---> try next char in S IF matchedChars == 0 THEN posInS := posInS + 1 ELSE matchedChars := F[matchedChars] // Seek back within pattern RETURN result This algorithm goes through the string once. With each character, it tries to match the beginning of P. When the start of P has been found, the algorithm tries to match the whole string. While doing so, counts the number of characters which already match. If the complete match fails, it gives up the complete match. But F tells it, how many characters it can nevertheless counted as found, because while scanning the incomplete pattern, a new pattern may have started in between. Invariant A says that counts the number of chars, which match from the beginning of P with the suffix of the currently read prefix of S. This is maintained trough all statements of the loop. Invariant B says that no occurences are left out. By invariant A, an occurence must be signalled by ==|P|. This case is cared for by the IF-statement. This algorithm runs in time O(|S|), because the difference of posInS - matchedChars increases with each main loop execution. Having precomputed F in time O(|P|), F can be used for a search in different strings. Now it becomes obvious, why F[i] considers the substring preceding i instead of the substring including position i: If, during the comparison, a character at position i mismatches, then one knows that the preceding substring was all right. If F[i] is used now, it must refer to this preceding substring. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Hashing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Key for a type t: A function key:t->tkey, which maps each instance of the type to an element of a key-set tkey, such that no two instances get the same key. List-Dictionary of type t: The type CLASS ListDictionary of T EXTENDS List of T key : FUNCTION T -> tkey M : Integer := 0 // Size PROCEDURE insert(e : T) this.pushFront(e) M := M + 1 PROCEDURE remove(k : tkey) IF this==nil THEN RETURN IF key(this.car)==k THEN this:=this.cdr M := M - 1 ELSE this.cdr.remove(k) FUNCTION find(k : tkey) : T IF this==nil RETURN nil IF key(this.car)==k THEN RETURN this.car RETURN this.cdr.find(k) Dictionary of type t: A type, which has operations equivalent to those of a list-dictionary of type t. Array-dictionary of type t: An array of type t, which has the operations of a dictionary. Hash-Function for an array-dictionary of type t: A function h:tkey->Natural, which maps an element of t given by its key to a position in the array. Collision conflict in an array-dictionary: The phenomenon that two elements are mapped to the same position by the hash-function. Closed Hash-Dictionary, Hash-Dictionary with Chaining of type t: An array dictionary of type t, which maintains a list for each array position. If two elements are mapped to the same position, they are simply stored both in the corresponding list. The definition for a Hash-Dictionary of type t with key key and hash-function h is the following: CLASS ChainingHash of type T m : Integer // Number of possible values of key : FUNCTION T -> tkey h : FUNCTION tkey -> {0,...m-1} d : array[this.m] of List of T M : Integer := 0 // Number of elements PROCEDURE insert(e : T) d[h(key(e))].pushFront(e) M := M + 1 PROCEDURE remove(k : tkey) d[h(key(e))].remove(k) // Using the List-Dictionary remove M := M - 1 FUNCTION d.find(k : tkey) : T RETURN d[h(key(e))].find(k) // Using the List-Dictionary find Inserting needs constant time, whereas removing and deleting may in the worst case take a time proportional to the number of elements in the list. But the expected time of removing and finding is O(1), if the size of the array d.size and the number of elements |d.M| in it are about the same. Proof: The expectation of the length of a list has to be determined. Be X:= d[h(k)].length for any key k, then X=|{e in d.M : h(key(e)) == h(k) }| We define the random variable X(e) := h(key(e))==h(k)?1:0 Then it holds E(X) = E(SUM {X(e) | e in d.M}) = SUM {E(X(e)) | e in d.M} = |d.M|*P(Xe=1) = |d.M|*1/d.size = |d.M|/d.size If the array has approximately the size of the number of elements, then the expected list length reduces to 1. // Warning: The notions "closed hashing" and "open hashing" are // sometimes used vice versa! Therefore, this summary uses the // non-ambiguous terms. Universal Hash-function family for an array-dictionary d: A set of hash-functions, where the probability of two elements being mapped to the same position is less than or equal to 1/d.size: P(h(key(e1))==h(key(e2))) =< 1/d.size "Vector-hash-function" for a dictionary d: The hash-function h[a](x[1],...x[k]) = SUM {x[i]*a[i] | i=0..k} mod d.size where * the key-function maps each element to a vector of {0,...d.size-1}^k. This is why the input to h is (x[1],...x[k]) * a=(a[1],...,a[k]) is a vector of {0,...d.size-1}^k, which characterizes h * d.size must be a prime number The set of all these hash-functions is universal. Proof: We consider two key-vectors x=(x[1],...x[k]) and y=(y[1],...y[k]), which differ in at least one position, say j: x[j]!=y[j]. There are {0,...d.size-1}^k different a and hence {0,...d.size-1}^k different h[a] in this family. We want to count those, which map x and y to the same position: h[a](x) = h[a](y) SUM {x[i]*a[i]} mod d.size = SUM {y[i]*a[i]} mod d.size 0 mod d.size = SUM {a[i]*(y[i]-x[i])} mod d.size a[j]*(y[j]-x[j]) mod d.size = SUM {a[i]*(y[i]-x[i]) | i=0..k, i!=j} mod d.size a[j] mod d.size = SUM {a[i]*(y[i]-x[i]) | i=0..k, i!=j} / (y[j]-x[j]) mod d.size Since a[j] {0,...2^20-1} h[t](x) := (x mod 2^20) xor t(x div 2^20) where t is any function t: {0,...2^12-1} -> {0,...2^20-1} Picture: We get a binary input x of 12 higher and 20 lower bits: xxx12xxxx xxxxxxx20xxxxxxxxxxx __ <- A := (x mod 2^20) | | xor '-- t --> ttttttt20ttttttttttt __| <- B := t(x div 2^12) | rrrrrrr20rrrrrrrrrrr <-' <- A xor B The set of all these hash-functions is universal. Proof: We have to show that for two different given keys x and y and a random t, the probability of x and y being mapped to the same cell is less than or equal to 1/d.size, where d.size=2^20: All x: All y: Random t: x!=y => P(h[t](x)==h[t](y)) There are two cases: * If the high parts of x and y are equal, then t(x div 2^20) and t(y div 2^20) will evaluate to the same result. Since the low parts of x and y are different, this difference will be passed directly to the result. Thus, the probability of x and y being mapped to the same cell is 0 in this case. * If the high parts are different, we need to find those t, which map x and y to the same cell, i.e. h[t](x) = h[t](y) <=> (x mod 2^20) xor t(x div 2^20) = (y mod 2^20) xor t(y div 2^20) We add an (y mod 2^20) left and right <=> (y mod 2^20) xor (x mod 2^20) xor t(x div 2^20) = (y mod 2^20) xor (y mod 2^20) xor t(y div 2^20) <=> (y mod 2^20) xor (x mod 2^20) xor t(x div 2^20) = t(y div 2^20) x and y are fixed. Hence, the left part of the equation is also fixed. Now x and y are mapped to the same cell, iff t(y div 2^20) takes exactly this fixed value. Since t is chosen randomly and t(y div 2^20) has the choice among 2^20 values, the probability is 1/2^20, which is 1/d.size. // Why is the left side of the equation fixed, as it contains t (?) Thus, the overall probability of the two keys being mapped to the same cell is less than 1/d.size, the family is universal. Open Hash-Dictionary, Hash-Dictionary with Probing: An array dictionary of type t, which seeks for a new free place in the array if an element is mapped to an occupied cell. Hash-Dictionary with Linear Probing: A hash-dictionary with probing, which uses the next free array cell whenever a new entry collides with a present entry. If the last array cell is bypassed, the search continues at the beginning. Hash-Dictionary with Bounded Linear Probing: A hash-dictionary with linear probing, which maintains additional space at the end of the array. Such a dictionary type t is defined as CLASS BoundedLinearProbingHash of T m : Integer // Number of possible values of key : FUNCTION T -> tkey h : FUNCTION tkey -> {0,...m-1} key : T size : Integer // Real size d : array[this.size] of T m' : Integer := size - d.m PROCEDURE insert(e : t) // The element should go here pos:=h(key(e)) // Seek for next free place WHILE this.d[pos]!=nil DO pos := pos+1 IF pos==this.size-1 THEN THROW Exception("Too many entries") // Put it in the free position this.d[pos]:=e FUNCTION find(k :tkey) : t // Should be here: pos:=h(k) // Seek pos and subsequent places WHILE this.d[pos]!=nil DO IF key(d[pos])==k THEN RETURN this.d[pos] pos := pos + 1 RETURN nil PROCEDURE remove(k : tkey) // First seek the element pos:=h(k) IF this.d[pos]==nil RETURN WHILE key(d[pos])!=k DO pos := pos + 1 IF this.d[pos]==nil THEN RETURN // We plan a hole for this.d[pos]. // Seek next element which could fill it pos' := pos + 1 WHILE this.d[pos']!=nil DO // Would this.d[pos'] have liked to go to this.d[pos]? IF h(key(d[pos'])) =< pos THEN // Put this.d[pos'] where it wants to be this.d[pos]:=d[pos'] pos := pos' pos' := pos' + 1 this.d[pos]:=nil Thus, there is never a hole between the actual position of an element and the would-like-to-be-position of an element: All j: j > h(key(d[i])) && j < i => d[j]!=nil ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Minimum Spanning trees ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Weight of a spanning forest of a weighted graph: The sum of the weights of its edges. WUC Graph: Weighted undirected and connected graph. Minimum spanning tree, MST of a WUC graph (V,E,c): A spanning tree, which has a minimal weight among all spanning trees of (V,E,c). c must be positive. Minimum spanning forest of a weighted undirected graph (V,E,c): A spanning forest, which has a minimal weight among all spanning forests of (V,E,c). Cut-Property: The fact that the lightest edge of a cut in a WUC graph can be used in a minimal spanning tree of that graph. Proof: Be (V,E,c) a WUC graph, be S a cut, i.e. a subset of the nodes. Be e* the minimal edge. Be R the connected component which is connected by e* to S. Either e* is in the MST, then the Cut-Property holds anyway. If e* is not in the MST, then any other edge e between R and S has to be in the MST. But since S and R are both connected components, both e and e* will do the job of connecting S and R. Hence, e can be replaced by e* in the MST and the Cut-Property holds. Cycle-Property: The fact that the heaviest edge on a cycle in a WUC graph is not needed in a MST. Proof: Suppose a MST uses the heaviest edge e* of a cycle. Then any other edge e of this cycle must miss in the MST. If e* is replaced by e in the MST, the result is still a MST. Jarnik-Prim-Algorithm, Prim-Algorithm: The following algorithm, which constructs a MST of a WUC graph: FUNCTION prim((V,E,c) : WUCGraph) : Set of Edge n := |V| // Edge to the predecessor of each node in the MST pred : array[n] of Edge // Distance of each node to the MST dist : array[n] of Integer := <+oo,...> // Nodes to be examined Q : Priority Queue of Node with priority dist := <> // Start with any node s0 := eps(V) Q.insert(s0) WHILE Q != <> DO s := Q.deleteMin() dist[s] := 0 FOREACH (s,v) in E DO IF c(s,v) < dist[v] THEN dist[v] := c(s,v) pred[v] := (s,v) IF v in Q THEN Q.decreaseKey(v) ELSE Q.insert(v) RETURN { pred[v] | v in V ,v != s0 } This algorithm works similarly to Dijkstra: It starts with a partial MST T, which initially only consists of an arbitrary start node s0. Then, it greedily adds the smalles edge leaving T to T. Thereby, T grows until all nodes have ben eaten. Similarly to Dijkstra, the algorithm uses time O(n+m) plus the time of the priority queue, which is O(m+n*log(n)) all in all with a Fibonaccy Priority Queue (not treated here) or O((m+n)*log(C)) with a Radix-Heap. Union-Find-Type: A type which can memorize a partition of a set of Integers M={1,...n}, where n is a given natural number. A partition is a set of subsets of M. In each subset, on element is chosen as a "leader", which represents the subset. For each element, one establishes parent-pointers, lead to the leader of the current subset. The type is CLASS UnionFind parent : array[n] of Integer := <1,2,3,...n> age : array[n] of Integer := <0,0,...> FUNCTION find(i : {1,...n}) : {1,...n} // Returns the leader of the subset, which i is in IF parent[i]==i THEN RETURN i i' := find(parent[i]) // Now as we found the leader, store it directly as a shortcut // in parent, avoiding the search for a future query // ("Path compression") parent[i] := i' RETURN i' PROCEDURE _link(i, j : {1,...n}) // Links two leaders to a new subset ASSERT i and j are leaders of different subsets // Have the "older" leader win IF age[i] follows for a given element the parent-pointer until it finds the subset-leader (the one which points to itself). In order to save this search for future calls, the parent pointer of the element is directly bent to the leader. The procedure connects two subsets to a new subset. Thereby, it just bends the pointer of the younger leader. The older leader then becomes even older. The age of a leader measures the size of its subset. By always chosing the older leader, it is guaranteed that those elements, which are in the greater (older) subset do not need an additional step in , whereas the elements in the smaller set do. But they are less, which improves the overall runtime. The runtime for all operations is nearly O(1). // The slides do not provide a proof for this runtime Union-Find-Datastructure: An instance of the Union-Find-Type. Kruskal-Algorithm: The following algorithm, which constructs a MST of a WUC graph: FUNCTION kruskal((V,E,c) : WUCGraph) : Set of Edges n := |V|; m := |E| T := UnionFind(n) result : Set of Edge := {} Esorted : array[|E|] of Edge := sortBy(E,c,INCREASING) FOR i:=0 TO |E|-1 DO (u,v) := Esorted[i] IF T.find(u) != T.find(v) THEN result := result \/ (u,v) T.union(u,v) RETURN result This algorithm grows the MST by joining subtrees. Each subtree is seen as a subset of V. These subsets are maintained in the Union-Find- Datastructure. First, every node is a tree of its own. The algorithm selects the lightest edge available and connects the two trivial subtrees to a new tree. It goes on like this, always chosing the next edge according to its weight: If it connects two subtrees, it is absorbed into the tree. If it is a connection within a tree, the cycle-property allows discarding the egde. Since and run in constant time, the dominating time factor is the sorting of the edges, which makes the overall runtime O(m*log(m)). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Geometry in the Plane ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // This is a section on two-dimensional objects. // All points and vectors are two-dimensional. // All sets of points defined here also form a type with the // appropriate components (e.g. Circle =(radius: Real, center: Point)). // All of these types have an operation "contains(p :Point)", // which determines whether p is part of the point set. // All real values or properties defined here in dependence upon // some other objects implicitely define a function // (e.g. "right turn: A sequence of Points x,y,z, such that..." // defines the function rightturn(x,y,z : Point) : Boolean). // These functions can be referred to in colloquial English // (e.g. "x, y and z form a right turn" instead of // rightturn(x,y,z)) Circular array, circular sequence: A sequence, the last element of which is its first element and all other element of which are distinct. For convenience, the last element is left away and one defines (x[1],...x[n]) := (x[2],...x[n],x[1]) x[i] := x[(n+i) mod n] Point: The type (x:Real,y:Real). This type inherits all operations of vectors. Coordinate of a point: One of its elements. Distance of two points p and q: The value |p-q|. Plane: The set of all two-dimensional points. One thinks of the plane as a blackboard, such that the notions "left" (small x-coordinate), "right" (large x-coordinate), "down" (small y-coordinate) and "up" (great y-coordinate) get meaning. Area: A subset of the plane. Epsilon-environment of a point p: The set of all points having a distance less than epsilon to p: { p' | |p-p'|0: Ex p': (p' not in area && |p-p'|0} is a subset of the area. Line segment between two points p and q: The set { p+r*(q-p) | 0 =< r =< 1} = { r*p+(1-r)*q | 0 =< r =< 1} Notation: pq with a bar above them, here: p--q It is p--q = q--p Arrow: A pair of points. Notation: pq with an arrow above them, here: p->q °: An abbreviation for "*pi/180" Angle of three points (p,q,r): The value angle(p,q,r) := (atn((r.y-q.y)/(r.x-q.x)) - atn((p.y-q.y)/(p.x-q.x)) + 2*pi) mod (2*pi) Example: p---------q 90° | r p---------q---r 180° r | p---------q 270° // This definition is quite arbitrary Turn matrix of three points p,q,r: The matrix 1 p.x p.y turn(p,q,r) := 1 q.x q.y 1 r.x r.y Left turn: A sequence of three points (p,q,r), such that angle(p,q,r)>180°. In this case, det(turn(p,q,r))>0. Right turn: A sequence of three points (p,q,r), such that angle(p,q,r)<180°. In this case, det(turn(p,q,r))<0. Collinear points: A sequence of three points (p,q,r), such that angle(p,q,r)==180° or angle(p,q,r)==0°. In this case, det(turn(p,q,r))==0. Point left of an arrow (p,q): A point r, such that angle(p,q,r)>180°. Point right of an arrow (p,q): A point r, such that angle(p,q,r)<180°. Intersection of two line segments: A point which belongs to both line segments. Given two line segments { a+r*(b-a) | 0 =< r =< 1} and { c+r*(d-c) | 0 =< r =< 1}, their intersection is the point a+r'*(b-a), such that 0=b and a and b lie on opposite sides of c->d. Line through the points p and q: The set { p+r*(q-p) | r in Real } = { r*p+(1-r)*q | r in Real } Line through the point p in the direction of vector d: The set { p+r*d | r in Real } This is the line through the points p and p+d. Line with the real parameters a, b, c: The set { (x,y) | a*x + b*y = c } This is the line through (0,c/b) and (1,(c-a)/b) and the line through (0,c/b) in the direction (1, -a/b ). Half-plane with the real parameters a,b,c: The set { (x,y) | a*x + b*y =< c } This is the set of points lying to one side of the line given by these parameters. Direction of a half-plane (a,b,c): The vector (-a, -b). It points to the direction where the half-plane is unbounded. Sequence of points in clockwise order: A sequence of points, each of which forms a right turn with its two successors. Sequence of points in counterclockwise order: A sequence of points, each of which forms a left turn with its two successors. Polygon vertices: A circular sequence P of points, such that no two line segments p[i]--p[i+1] p[j]--p[j+1] intersect (except for equal points) and such that no point is collinear with its two sucessors. Edge of polygon vertices: A line segment of a vertex and its successor. Polygon of polygon vertices P: The set { q | | (+oo,+oo)--q /\ E | mod 2 == 1 }, where E is the set of all edge points of P. That is: If I take a point exterior to the boundary and draw a line to q, then this line crosses the boundary an odd number of times. A polygon is mostly given by its vertices. We define the type Polygon as being a circular array. Circle with radius r around a center point p: A set of points, all of which have the distance r to p: { q | |q-p|==r } A circle is mostly given by its radius and center. Disc of a circle C: The set { q | | (+oo,+oo)--q /\ C | mod 2 == 1 }, that is the "inner area" of the circle. A disc is mostly given by its circle. Interior of a circle: Its disc. Convex area: An area S, such that for any two points p and q in S, the line segment p--q is a subset of S. This corresponds to an area without "inward corners". The intersection of convex areas is again convex. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Convex Hull ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Site: A marked point. This notion is used to distinguish certain given points from the other infinitly many points of the plane. Convex hull of a set of sites S: The smallest convex area containing S. If one thinks of the sites as nails in a board, then the convex hull corresponds to a rubber band (Gummiband) embracing all of them. If S is finite, then the convex hull is a convex polygon, which is mostly given by its vertices in clockwise order. Notation: CH(S) Reduction of the sorting problem to the convex hull problem: The following algorithm, which solves a sorting problem if a convex hull problem can be solved: FUNCTION sortByCH(A : array of Real) : array of Real points : array[A.length] of Point FOR i:=0 TO A.length-1 DO points[i] := (A[i],-A[i]^2) ch : circular array[A.length] of Point := convexHull(points) minx := ch[0].x FOR i:=0 TO A.length-1 DO IF ch[i].xq THEN counter := counter + 1 IF counter==|P|-2 THEN edges := edges \/ (p,q) result := array[edges.length] of Point result[0]:= eps(edges).startNode FOR i:=1 TO |edges|-1 DO result[i] := eps({ q | (result[i-1],q) in edges}) RETURN result This algorithm considers each pair of sites and checks whether all other sites lie right of the line egment defined by the pair. If so, it considers the pair an edge of the convex hull. Last, the set of edges is transformed to the sequence of vertices. This algorithm runs in time O(|P|^3). Incremental algorithm: An algorithm, which solves aproblem first for one constraint and then adds the other constraints one after the other. Upper Hull of a set of sites: The sub-sequence of the convex hull polygon vertices, in which each site has a larger x-coordinate that its predecessor. Graham's Scan: The following incremental algorithm, which computes the convex hull of a set of sites P: FUNCTION halfHull(P : array of Point, X : {left, right}) : array of Point result : List of Point := <> result.push(P[0]) result.push(P[1]) FOR i:=2 TO P.length-1 DO WHILE result.get(result.length-2), result.get(result.length-1), P[i] do not form a X turn DO result.popBack() result.pushBack(P[i]) RETURN Array(result) FUNCTION graham(P : Set of Point) : Polygon P' : array[|P|] of Point := sortBy(P,x,INCREASING) upperHull := halfHull(P', right) lowerHull := halfHull(P', left) RETURN Polygon(upperHull + reverse(lowerHull)) This algorithm first sorts the sites from left to right, i.e. by their x-coordinate. It considers the boundary of the convex hull polygon of P to consist of two parts: The "upper" hull (only right turns from left to right) and the "lower hull" (only left turns from left to right). Both are calculated by adding one point after the other to the hull half. If the new site does not form the correct turn with its two predecessors, then the predecessor is kicked out. Then the two hulls are connected. This algorithm runs in time O(n*log(n)). b_______c e Adding e will remove f of the / \ / convex hull a/ \/ f Jarvi's March: The following algorithm, which computes the convex hull of a set of sites P: FUNCTION jarvi(P : Set of Point) : Polygon result : List of Point := <> p0 : Point := (0,+oo) // Pick the lowest point p1 : Point := lowestPoint(P) P := P \ lowestPoint(P) // Loop until we are back at the beginning REPEAT // Find p2, such that angle(p0,p1,p2) is maximal p2 := eps(P) FOREACH p in P DO IF angle(p0,p1,p) > angle(p0,p1,p2) THEN p2 := p P := P \ {p2} (p0,p1) := (p1,p2) UNTIL p1==result.get(0) RETURN Polygon(result) This algorithm first picks the lowest site p1 of the set. Then it checks which other site makes up a maximal angle with p1--(0,+oo). This site is chosen as a next site in the vertex sequence. Then it checks again which remaining site makes up a maximal angle with the newly added site and the predecessor of it. This continues until the start site point has been met. If n is the number of sites and h is the number of convex-hull-vertices, then this algorithm runs in time O(n*h). Tangent of a point p to a connected area A: A line through p, such that the line has exactly one point q in common with A and the whole area lies to the right (or left) of p->q. The left tangent of a convex polygon P (i.e. where the polygon will lie to the right of it) can be calculated in time O(log(P.length)) as follows: FUNCTION leftTangentTouchPoint(P : Polygon, p1:Point) : Integer ASSERT p1 lies to the left of the polygon, P convex // phi(i) is the angle which P[i] forms with p1 and a // point in the negative y-infinity phi : FUNCTION Integer -> Real, phi(i) := angle(P[i],p,(0,-oo)) // We need to find the maximum of phi // Define some useful functions on phi phiGrows : FUNCTION Integer -> Boolean, phiGrows(i) := phi(i-1) < phi(i) < phi(i+1) phiShrinks : FUNCTION Integer -> Boolean, phiShrinks(i) := phi(i-1) > phi(i) > phi(i+1) isMax : FUNCTION Integer -> Boolean, isMax(i) := phi(i-1) < phi(i) > phi(i+1) isMin : FUNCTION Integer -> Boolean, isMin(i) := phi(i-1) > phi(i) < phi(i+1) // When does P[i]...P[j] contain the maximum? containsMax : FUNCTION Integer x Integer -> Boolean containsMax(i,j) := phiGrows(i) && phiShrinks(j) || phiGrows(i) && phiGrows(j) && phi(i)>phi(j) || phiShrinks(i) && phiShrinks(j) && phi(i) i (circular) We look at two indices i and j. By the help of phi(i) and phi(j), we can tell whether the maximum lies between i and j: If phi(i) decreases and phi(j) increases, then the maximum must be in between. The maximum is also in between, if both phi(i) and phi(j) increase, but phi(i) > phi(j). A number of other cases applies similarly. Now we do a binary search with i and j: If the maximum is in between, we look at the index k:=(i+j)/2. If the maximum is not in between, we swap i and j. Since the difference i-j is always divided by two, the total running time is O(P.length). // This is David Steurer's solution. Thanx! Chan's algorithm: The following algorithm, which computes the convex hull of a set of sites P: FUNCTION partialHull(P : Set of Point, m : Integer) : Polygon r := |P| / m // Partition P into r sets of at most m elements PP : array[r] of Set of Point := <{},...{}> FOR i:=0 TO r-1 DO FOR j:=0 TO m DO IF P=={} BREAK PP[i] := PP[i] \/ {pickAny(P)} // Compute the convex hulls of each subset PPhulls : array[r][m] of Point := // Do Jarvi's March on polygons p0 := (0,+oo) p1 := lowestPoint(P) p2 := p0 result : List of Point := FOR i:=1 TO m DO // Find p2, such that angle(p0,p1,p2) is maximal FOR j:=1 TO r DO left := PPhulls[j][leftTangentTouchPoint(p1,PPhulls[j])] right :=PPhulls[j][rightTangentTouchPoint(p1,PPhulls[j])] IF angle(p0,p1,left) > angle(p0,p1,p2) THEN p2 := left IF angle(p0,p1,right) > angle(p0,p1,p2) THEN p2 := right IF p2==lowestPoint(P) THEN RETURN Polygon(result) result.pushBack(p2) (p0,p1) := (p1,p2) RETURN nil FUNCTION chan(P : Set of Point) : Polygon t := 1 REPEAT m := min(2^(2^t), |P|) result := partialHull(P,m) UNTIL result != nil RETURN result The partialHull function partitions the set of sites into subsets of size m and computes the convex hull of each of them. By Grahams scan, this takes time n/m * O(m*log(m)) = O(n*log(m)). Then, it runs Jarvi's March on these polygons. Given that leftTangentPoint runs in time O(log(m)), this takes time O(h*n/m*log(m)), where h is the number of vertices of the final convex hull. If m=h, then the algorithm runs in time O(n*log(h)). Since h is unknown, it is found by repeated squaring in the chan-function: First m=4 is tried, then m=16, m=256, m=65536 etc. The t-th iteration takes time O(n*log(2^(2^t))) = O(n*2^t). If m>h, then the partialHull-function will find the hull, so the iteration is executed at most log(log(t)) times. The overall time amounts to SUM {n*2^t | t=1...log(log(h))} in O(n*log(h)). _ _ |:| Compute the small polygons by Graham, then /:\ |_| _ compute an all-enclosing polygon by Jarvi \:/ /:\ \:/ Polygon-merging: The following algorithm, which merges two disjoint convex polygons to a new one: FUNCTION polygonMerge(P1, P2 : Polygon) : Polygon ASSERT P1 and P2 convex and disjoint, P1 left of P2 // Find upper tangent i1 : Integer := rightMostPointIndex(P1) i2 : Integer := leftMostPointIndex(P2) WHILE not (P2 is right of P1[i1]->P2[i2] && P1 is right of P1[i1]->P2[i2]) DO WHILE not P2 is right of P1[i1]->P2[i2] DO i2 := i2 + 1 WHILE not P1 is right of P1[i1]->P2[i2] DO i1 := i1 - 1 // Find lower tangent j1 : Integer := rightMostPointIndex(P1) j2 : Integer := leftMostPointIndex(P2) WHILE not (P2 is left of P1[j1]->P2[j2] && P1 is left of P1[j1]->P2[j2]) DO WHILE not P2 is left of P1[j1]->P2[j2] DO j2 := j2 + 1 WHILE not P1 is left of P1[j1]->P2[j2] DO j1 := j1 - 1 // Merge the halves and tangents to the result result : List of Point REPEAT result.pushBack(P[j1]) j1 := j1 + 1 UNTIL j1==i1 REPEAT result.pushBack(P[i2]) i2 := i2 + 1 UNTIL i2==j2 RETURN Polygon(result) This algorithm assumes P1 to be left of P2. We find the rightmost vertex i1 of P1 and the leftmost vertex i2 of P2. They so to speak face each other. Then we move i1 and i2 up alternatedly, until we found the upper tangent. The condition "P1 lies right of P[i1]->P[i2]" can be checked by looking if both P1[i1-1] and P1[i1+1] lie to the right of P[i1]->P[i2]. We do the same to find the lower tangent and finally merge the left half of P1, the upper tangent, the right half of P2 and the lower tangent. The runtime is in O(n), where n is the total number of vertices. _______c z________ The run will be a--x, a--y, a--z, / \b y/ / b--z, c--z, done | | | / \ /a x\_____/ \_____/ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Line segment intersection ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Sweep line: An arrow (i,-oo)->(i,+oo), which is moving across the plane with increasing i. Sweep line techniques are used for solving problems in the plane. Thereby, the problem has already been solved to the left of the sweep line. Event point of a sweep line: A site, which will be part of the sweep line at some time. The sweep line always jumps from one event point to the next. Dead area of a sweep line: An area left of it. Active area of a sweep line: An area which has points in common with the sweep line. Dormant area of a sweep line: An area to the right of the sweep line. Line intersection sweeping: The following algorithm for calculating the set of all intersection points of a set of line segments: FUNCTION lineIntersect(L : Set of LineSegment) : Set of Point result : Set of Point := {} CONST INTERSECT : Point := (+oo,+oo) xstructure : Heap of (point:Point,to:Point) with priority point.x // Insert all end and start points of all line segments // into the xstructure FOREACH p--q in L DO xstructure.insert(( p ,q)) xstructure.insert(( q ,p)) // All line segments intersecting the sweep line ystructure : Heap of (point:Point,to:Point) with priority -point.y WHILE !xstructure.empty() DO (p,q) := xstructure.deleteMin() pred := ystructure.pred(p,q) succ := ystructure.succ(p,q) // We have an intersection point IF q==INTERSECT THEN result := result \/ {p} ystructure.delete((p,q)) ppred := ystructure.pred(pred) ssucc := ystructure.succ(succ) IF ppred intersects with succ to the right of the sweep line (p.x,-oo)->(p.x,+oo) THEN xstructure.insert((intersectPoint(ppred,succ),INTERSECT)) IF ssucc intersects with pred to the right of the sweep line (p.x,-oo)->(p.x,+oo) THEN xstructure.insert((intersectPoint(ssucc,pred),INTERSECT)) ystructure.swap(pred,succ) CONTINUE WHILE // We have a line start point IF p.xq.x THEN IF pred intersects with succ to the right of the sweep line (p.x,-oo)->(p.x,+oo) THEN xstructure.insert((intersectPoint(pred,succ),INTERSECT)) ystructure.delete(p,q) RETURN result Assume no line segment to be vertical and no two line endpoints or intersection points have the same x-coordinate. This algorithm moves a sweep line across the plane. All line segment start and end points (left and right points) are event points. These event points are stored left to right in a heap called xstructure. Furthermore, all intersection points are also event points. Since these are unknown at the beginning, they are calculated on the fly as the sweep line moves and added to xstructure. All time, a heap called ystructure stores all line segments intersecting the sweep line, sorted from above to below. Now the sweep line moves to the next event point and different things happen: * If we meet the start of a line, then this line is added to ystructure and eventual intersections of this line with its neighbors in ystructure are recorded as event points in xstructure * If we meet the end of a line, then we check whether its neighbors in the ystructure intersect to the right of the sweep line. If this is the case, we add the intersection point as a new event point to the xstructure. Then we delete the line from the ystructure. | i\| When line j ends, i and k are checked | / for interection and q is found and |\ / added to the x-structure j_____| \/.....q | /\ |/ \ | \ k/| | sweep line * If we met an intersection point, we check whether the upper line will intersect with the successor of the lower line in the ystructure. We also check whether the lower line will intersect with the predecessor of the upper line in the ystructure. We add the intersection points as event points to the xstructure. Then we interchange the upper and the lower line in the ystructure. | i\ | / When the intersection of i and k is met, j_____\__|__/______ i and k swap their order in the \ | / y-structure. It is checked whether the \|/ lower line k meets the predecessor of the X upper line i, which is j here /|\ / | \ k/ | |sweep line Be n the number of lines and s the number of intersection points. As each of the heap operations takes time O(log(n)), the total time amounts to O((n+s)*log(n)). Proof of correctness: Be x an intersection point of two line segments p--q, r--s. Assume that p is left of q and r is left of s. All points p,q,r,s will be event points. Assume p is left of r, so p will be in the ystructure when r is inserted. Assume all intersection points left of the sweep line have already been found (trivially true in the beginning). There are three cases: * There is no line segment reaching into the polygon p,x,r. As there is no line in between, p and r will be neighbors in the ystructure. When r is inserted, the intersection with its neighbors is checked and so x will be found. * There are line segments reaching into p,x,r. Be k the end point of the rightmost of these lines. As soon as k is deleted from the xstructure, it is also deleted from the ystructure. p and r become neighbors and their intersection is checked, x is found. * There is a line segment intersecting with p--q left of x and to the right of the rightmost endpoint of any line reaching into p,x,r. Then this intersection point will have been found and will be an event point. Before exchanging the two lines in the ystructure, eventual intersections with the neighbors are checked and x is found. If x is found, it is added as an event point to the xstructure and reported as an intersection. Polygon intersection sweeping: A modified version of the line intersection algorithm, which calculates the intersection of two convex polygons. The edges of the polygons are the line segments the algorithm works on. The y-structure can only contain 4 line segements at a maximum. The sweep line starts at the leftmost point of both polygons. The xstructure does not contain all line end- and start-points, but just the end-points of the line segments stored in the ystructure. Then there is also a fixed number of possible intersection points. As a result, the xstructure also contains a constant number of lines. Every endpoint is at the same time a startpoint. New endpoints are added on the fly. Hence all heap accesses take constant time. Be n the number of vertices and s the number of intersections of the polygon edges, then the total time is O((n+s)*1) = O(n). Knowing the intersection points in the order from left to right, the polygon intersection can be calculated. Half-Plane-Intersection-Algorithm: The following algorithm, which, given a set of half-planes, determines their intersection: FUNCTION hpintersect(H : Array of HalfPlane) : Set of LineSegment IF |H|==1 THEN RETURN { LineSegment(Line(H[0])) } n2 := H.length/2 RETURN polygonIntersection_(hpintersect(H[0..n2]), hpintersect(H[n2..H.length-1]) This is a recursive divide-and-conquer-approach, which relies on the fact that the polygon-intersection sweeping algorithm also works for "unbounded polynoms", if adjusted appropriately. If n is the number of half-planes, then we have log(n) recursive calls, each of which takes time O(n) for the intersection, resulting in a total time of O(n*log(n)). Binary tree: A tree, all nodes of which have 0, 1 or 2 children. Check-tree: The following type CLASS CheckTree n : Integer depth : Integer tree : array[depth][n] of (mean:Real,checkedLeaves:Integer) PROCEDURE check(r: Real, b:Boolean) // Checks a leaf pos:=0 FOR level:=0 TO depth-1 DO tree[level][pos].checkedLeaves +:= b?1:-1 IF r>tree[level][pos].mean THEN pos:=pos*2+1 ELSE pos:=pos*2 FUNCTION inbetween(r,s:Real):Integer // Returns the number of checked leaves betwen r and s, inclusive result:=0 pos:=0 FOR level:=0 TO depth-2 DO IF r>tree[level][pos].mean THEN // Go right pos:=pos*2+1 ELSE // Go left, add right checked leaves pos:=pos*2 result:=result+tree[level+1][pos*2+1].checkedLeaves IF tree[depth-1][pos].mean==r THEN result := result + tree[depth-1][pos].checkedLeaves FOR level:=0 TO depth-2 DO IF s>tree[level][pos].mean THEN pos:=pos*2+1 ELSE pos:=pos*2 result:=result-tree[level+1][pos*2+1].checkedLeaves RETURN result FUNCTION makeCheckTree(S : Set of Real) // Returns a CheckTree for S t : FakeTree t.depth:= ceil(log2(|S|)) t.n:=2^t.depth // Sort the reals sorted : array[t.n] of Real := sort(S,INCREASING) // The sorted numbers are the leaves FOR i:=0 TO t.n-1 DO result.tree[t.depth-1][i]:=(sorted[i],0) // Create a tree over them FOR level := t.depth-2 TO 0 DO FOR node := 0 TO 2^level DO t.tree[level][node].mean := (t.tree[level+1][node*2]+ t.tree[level+1][node*2+1])/2 t.tree[level][node].checkedLeaves := 0 RETURN t This CheckTree-class serves only for the following Orthogonal-Line- Sweeping-algorithm. It can handle a set of reals S, some of which are "checked". It sorts all numbers of S and builds a binary tree over them. Each node of the binary tree stores the mean value of its children, so that a leaf can be found by passing from the root to the leaves. Each node also stores the number of checked leaves below it. If a number is to be checked or unchecked, the corresponding leaf is found and the number of checked leaves is decreased on the way from the root to the leaf. If the number of checked reals between two reals r and s is to be determined, we run down the tree from the root to r and sum up all checked leaves of all right children, whenever we chose a left child. Thereby, we count all checked leaves right to r. Now we do the same again for s and sum up all checked leaves right to s. The difference is the number of intermediate checked leaves. All operations take time O(log(|S|)). Orthogonal-Line-Sweeping: The following algorithm, which calculates the number of intersections between a set of horizontal line segments and a set of vertical line segments: FUNCTION OlineIntersect(H,V : Set of LineSegment) : Integer result : Integer := 0 xstructure : PriorityQueue of (from, to:Point) with priority from.x // Insert all line points into the xstructure FOREACH p--q in H \/ V DO xstructure.insert((p,q)) xstructure.insert((q,p)) // All horizontal line segments will be checked ystructure : makeCheckTree({y | y is a y-coordinate in H}) WHILE !xstructure.empty() DO (p,q) := xstructure.deleteMin() // H-line start IF p.xq.x THEN ystructure.check(p.y,false) // V-line IF p.x==q.x THEN result := result + ystructure.inbetween(p.y,q.y) RETURN result This algorithm moves a sweep line across the plane like the line intersection sweeping algorithm. All line points are event points. When a horizontal line is met, its y-coordinate is stored in the ystructure. This entry is deleted, when the horizontal line is over. When a vertical line is met, the number of horizontal lines between start- and end-point of the vertical line is added up. The time is O((|V|+|H|)*log(|V|+|H|)). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Linear Programs ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ LP, Linear Program: A pair of * A maximize expression: c*x * A set of constraints d[1]*x =< d'[1] d[2]*x =< d'[2] ... where x is a vector variable, the d[i] are constant vectors and the d'[i] are real constants. The dimension of x is called the dimension of the linear program. Here, we will deal with two-dimensional linear programs. The following variations are allowed: * One may also use the relation '>=', since it can be expressed by multiplying the line by -1 and using '=<'. * One may also use the relation '=', as it can be expressed by a line with '=<' and a line with '>='. * One may express the lines as sums of the form SUM {d[i][j]*x[j] | i=1..} =< d'[i] Feasible region of a linear program: The set of all points (x[0],x[1]), which satisfy the constraints of the linear program. Pointer of a linear program: The vector (c[1],c[2]). It holds: For any point p, p+v has a greater value for the maximize expression than p. The value of a point (x,y) for the maximize expression is c[1]*x + c[2]*y = (x,y)*(c[1],c[2]) Solution of a linear program: The point, which lies in the feasible region and the coordinates of which maximize the maximize expression. If the feasible region is empty, a solution does not exist. If the feasible region is unbounded in the direction of the pointer, the solution is infitity. A solution of a linear program can be found in polynomial time. A typical execution time is O((m+n)*N), where n is the dimension of the linear program, m is the number of constraints and N is the number of non-zero components of the vectors d[i]. // Algorithm treated later Line-adding-algorithm: The following incremental algorithm, which solves a linear program: FUNCTION linProg(H : Array[n] of HalfPlane, pointer : Vector) : Point // Find two halfplanes, which bound the feasible region in the // direction of p, put them to the beginning of H j := 0 FOR i:=0 TO n-1 DO // If there is no upper bound, we're done IF i==H.length THEN RETURN +oo*pointer Be ph a point on the boundary of H[i] IF not ph+pointer in H[i] THEN // h bounds feasible region swap(H[j],H[i]) j := j + 1 IF j==2 THEN BREAK FOREACH // tentative result: Intersection of these two upper bounds result : Point := intersection(Line(H[0]),Line(H[1])) // Store for every half-plane its boundary with the feasible // region by two points and . We order these points, such // that a*pointer > b*pointer. Hboundary : array[n] of (a,b:Point) Hboundary[0]:= (result, ) Hboundary[1]:= (result, ) // Add other half-planes, check whether result is not OK FOR i:=2 TO n-1 DO IF not result in H[i] THEN // Check intersection with H[i] with all previous ones bestsofar := -oo Hboundary[i]:=(, ) FOR j:=0 TO i-1 DO intsct : Point = intersection(Line(H[i]),Line(H[j])) IF intsct in Hboundary[j].a--Hboundary[j].b && intsct*pointer < result*pointer THEN bestsofar := intsct*pointer result := intsct Hboundary[j]:= (intsct,Hboundary[j].b) IF intsct*pointer > Hboundary[i].a*pointer THEN Hboundary[i]:=(intsct,Hboundary[i].b) ELSE Hboundary[i]:=(Hboundary[i].a,intsct) // If we didn't find anything, there is no solution IF not result in H[i] THEN RETURN <> RETURN result This algorithm first finds two half-planes, which bound the feasible region in direction of the pointer. If there is none, the problem is solved. We take the intersection point of the upper bounds as a tentative result. Then we add the other half-planes one after the other. If the tentative result point lies within the newly added half-plane, we can keep the point. Else, we check all intersections of the new line with the previous lines, as these may become new results. We pick the best among them. We only take those intersection points into account, which lie on a boundary with the feasible region. In the best case, this algorithm has a run-time of O(n): The lines are added in such an order, that the seek for intersection points is never necessary. In the worst-case, it has run-time of O(n^2), because each of the newly added lines may require a check with all previous lines. i\ /j For {i,j}, p was the optimum. Now k is \/......p added. Since p does not lie in the half- ^ /\ space blow k, the intersections of k with | c ___/__\_____k all previous lines have to be checked. / \ Shuffling: The following algorithm: FUNCTION shuffle(A: array) : array FOR i := A.length-1 DOWNTO 1 DO swap(A[i],A[random(i)]) // random(i) in {0...i-1} RETURN A This algorithm creates a randomized permutation of the array. The probability of any given element being at any given position is 1/A.length. Randomized algorithm: An algorithm, which, for the same inputs, may do different things. Expected runtime of a randomized algorithm: Its average run-time for an infinite number of runs on the same inputs. Backward-analysis: The following technique for determining the run-time of an incremental randomized algorithm: Instead of adding a constraint and asking what happens with what probability, one assumes a situation with the first constraints already given and asks what happens with what probability if the last constraint is removed. Randomized line-adding-algorithm: The line-algorithm, which shuffles the array of half-planes first. This shuffling excludes the first two half-plains, which have been found to be upper bounds. The shuffling causes the algorithm to run in expected time O(n): If the half-plains are added in random order, chances are good that * either the tentative result point already satisfies the newly added constraint or * if it does not, there are few previous lines with which to check intersections. Proof: We do a backward analysis. Assume the half-planes H[0]...H[i] to be given and p to be the tentative result point. Now we remove H[i] and imagine the situation before H[i] was added. Was the previous p different from the new p? p always lies on an intersection of two lines. These lines are said to determine p. If by chosing H[i], we incidentally chose a line which determined p, then this means p was different before H[i] was added. What is the probability that by removing H[i], we remove a line determining p? Only two lines out of the i lines determine p, hence the probability is 2/(i-2). (We subtract 2 to take into account the first two predetermined lines). With this probability of 2/(i-2), we will have to check all (i-1) previous lines for intersection points with H[i]. The expected running time is thus O(n)+ for finding the two upper bound lines O(SUM {2/(i-2)*(i-1) | i=3..n}) for the intersection checks = O(n). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Smallest enclosing Disc & Duality ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Smallest enclosing disc of a set of sites: A disc, which contains all sites and has the smallest radius of all these discs. This disc can be computed by the following incremental algorithm: FUNCTION smallestDisc3(q1, q2, q3 : Point) : Disc // Computes the smallest disc containing q1, q2, q3 p : Point := (q1+q2+q3)/3 RETURN (radius: |q1-p|, center: p) FUNCTION smallestDisc2(P : array of Point, q1,q2 : Point) : Disc // Computes the smallest disc containing P and with // q1, q2 on its boundary D : Disc := (radius:|q1-q2|, center:(q1+q2)/2) FOR i := 0 TO P.length-1 IF P[i] is not contained in D THEN D := smallestDisc3(q1,q2,P[i]) RETURN D FUNCTION smallestDisc1(P : array of Point, q1 : Point) : Disc // Computes the smallest disc containing P and with // q1 on its boundary D : Disc := (radius:|q1-P[0]|, center:(q1+P[0])/2) FOR i := 1 TO P.length-1 IF P[i] is not contained in D THEN D := smallestDisc2(P[0..i-1],P[i]) RETURN D FUNCTION smallestDisc(P : array of Point) : Disc P := shuffle(P) // Preliminary result D : Disc := (radius:|P[0]-P[1]|/2, center:(P[0]+P[1])/2) FOR i:=2 TO P.length-1 DO IF P[i] is not contained in D THEN // D is the disc containing P[0..i-1] with P[i] // on its boundary D := smallestDisc1(P[0..i-1], P[i]) RETURN D This algorithm starts off with a function to compute the smallest disc of three sites. The next function computes the smallest disc for a set of sites, where it is already known that two sites q1 and q2 lie on its boundary. It just computes the smallest disc of q1,q2 and the first site of the set. If any other site p is not in this disc, then the disc around q1,q2,p is taken instead. Does this work(?) Suppose q1 p1 p2 q2 First, the circle goes through q1 and q2. It contains p1, so the tentative circle does not change for p1. Now, p2 is added and the circle is computed as smallestDisc3(q1,q2,p2), which may grow so large, that the circle line between q1 and q2 might become straight, p1 is excluded from the disc! This doubts lemma 1(iii). The next function computes a smallest enclosing disc for a set of sites, where it is already known that one site q1 lies on its boundary. The sites of the set are added incrementally and whichever site is not in the disc must be on its boundary. In this case, the previous function is called. The main function does a shuffling and proceeds according to the same pattern. The run-times are as follows: * smallestDisc3: O(1) * smallestDisc2 with i points: O(i) * smallestDisc1 with i points: It has to be calculated how often smallestDisc1 calls smallestDisc2. smallestDisc2 is called for j-1 points, if P[j] does not lie on the tentative disc. We do a backwards analysis: We suppose a given state j with a disc containing all sites. Now we remove the last site P[j]. What is the probability that the disc was different before? The smallest disc was only different if P[j] is one of the three sites on the boarder of the current disc. This is the case with a probability of 3/j. Since in this case, smallestDisc2 will have been called before, we can compute the running time of smallestDisc1 as SUM {O(j)*3/j + 1*(1-3/j) | j=2..i } = O(i) * smallestDisc with n points: O(n) by a similar analysis Dual point of a line y=a*x-b: The point(a,b). Notation for the dual point of a line L: L* Dual line of a point (a,b): The line y=a*x-b. Notation for the dual line of a point P: P* Self-inverse-lemma: The fact that (p*)=p. Proof: (p*)* = ({(x,y)|y=p.x*x-p.y})* = (px.py) = p Order-reversing-lemma: The fact that iff point p lies above line L, then point L* lies also above line p*. This also applies analogously to the inverted situation and if p lies on L. Proof: Suppose p=(p.x,p.y) lies above L:y=la*x-lb p.y>la*p.x-lb <=> lb > p.x*la-p.y ...which means that L* lies above p* Incidence-preserving-lemma: The fact that the lines L1 and L2 intersect at point p, iff the line p* passes through the points L1* and L2*. Proof: Suppose lines L1 and L2 intersect in point p. p is on L1 && p is on L2 <=> L1* is on p* && L2* is on p* (by order-reversing-lemma) Thus p* passes through L1* and L2*. Collinearity/Coincidence-lemma: The fact that three points are collinear, iff their dual lines intersect in a common point. If the three points are exactly one above the other, the lemma only holds if one accepts parallel lines to intersect in (+oo,+oo). Proof: Suppose lines L1, L2, L3 intersect in point p. p is on L1 && p is on L2 && p is on L3 <=> L1* is on p* && L2* i on p* && L3* is on p* (by order-reversing-lemma) Thus, L1*, L2* and L3* are collinear. Upper envelope of lines: The boundary of intersection of the half-planes above these lines. LowerHull/UpperEnvelope-Theorem: The fact that the left-to-right order of the lower hull of a set of sites is the same as the left-to-right-order of the sequence of their dual lines on their upper envelope. Proof: Suppose p[i]--p[j] is in the lower hull of a set P of sites. all p lie above p[i]--p[j] <=> all p* lie below (p[i]--p[j])* (by order-reversing-lemma) <=> all p* lie below p[i]* /\ p[j]* (by incidence-preserving-lemma) <=> all dual lines below intersection of p[i]* and p[j]* This means that p[i]* and p[j]* form two adjacent edges of the upper envelope. Suppose p[i] is left of p[j]. Then p[i]* will have a smaller slope than p[j]*. Hence the contribution of p[i]* to the upper envelope must lie left of the contribution of p[j]*. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Voronoi-Diagrams ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Binary tree structure of type T with priority f:T->Real: A heap of type T with priority f:T->Real with all its operations, basing on the following type structure: CLASS BinaryTree of type T with priority f:T->Real left, right : BinaryTree of type T with priority f:T->Real element : T Balanced binary tree: A binary tree, where no leaf is on a level smaller than the maximum level - 2 (?). Balanced binary tree structure: A binary tree structure, which with each delete- or insert-operation re-organizes the tree such that the corresponding binary tree is balanced. As a result, all operations always take time O(log(n)), if n is the number of nodes in the tree. Box with the edge points p and q: The set p--(p.x,q.y) \/ p--(p.y,q.x) \/ q--(p.x,q.y) \/ q--(p.y,q.x). Non-crossing line segments: Two line segments, which either have no point in common or just one end-point. Curve: A sequence of (possibly infinitely small) non-crossing line segments, each of which has the start-point identical to the end-point of its predecessor -- except for the first one. Planar graph: An undirected graph, the nodes of which are points and the links of which are curves between the points. The curves may only have the nodes in common. Straight Planar graph: A planar graph, all edges of which are line segments. There are additional line-segments (edges) allowed, which start at one node and are unbound in one direction. These lines are usually dealt with by adding a big box around the graph, such that these lines can be given by their start point and their intersection point with the box. Face of a straight planar graph: A (possibly unbounded) polygon of its nodes, which does not include other nodes. Euler's formula: The fact that in a straight planar graph the number of nodes plus the number of faces minus the number of edges equals two: nodes+faces-edges=2 This implies that the expected degree of any vertex is 6. (?) Cocircular points: Points, for which there is a circle, such that they all lie on it. General position assumption for a set of points P: The assumption that no 4 points of P are cocircular. Voronoi-cell of a site q in a set of sites P: The set of all points in the plane, which are closer to q than to another site. V(q) := {x | |q-x| =< |p[i]-x| for all p[i] in P} A voronoi-cell is a possibly unbounded convex polygon and given by the sequence of its vertices. Voronoi-diagram of a set P of sites: The straight planar graph made up by the set of all Voronoi-cells of P. The general position assumption will be made. The following properties hold: * Each point on an edge between two Voronoi-cells V(p[i]) and V(P[j]) has the same distance to the sites p[i] and p[j] * A point, which has the same distance to three sites is a vertex of the Voronoi-diagram. It is the center of the circle passing through the three sites with no other site inside the circle. * The vertices of the Voronoi-diagram all have degree 3. * The Voronoi-diagram has |P|=:n faces, less than 2*n-5 vertices and less than 3*n-6 edges. Notation: Vor(P) \ | . / \ . |__/ . \___/ . \_______ . / . \___/ . \___/ . \ / . | \ / | Naive Voronoi algorithm: An algorithm, which computes the Voronoi-diagram of n sites in time O(n^2*log(n)) by computing the intersection of n half-planes per site. (?) Distance of a site to a line: The minimum among the distances of all points on the line to the site. Parabola of a site s and a sweep-line: The set of all points of the plane, which have the same distance to s as to the sweep-line. If we have a horizontal sweep-line at y-coordinate ly, then the parabola is given by {(x,y) | (y-ly)^2 = (x-s.x)^2+(y-s.y)^2} It is assumed that the sweep line moves downwards. Beach line of sites S and a horizontal sweep-line: The set of all those points on the parabolae, which do not have another parabola point below them. Breakpoint of a beach line: A point of a beach line, which lies on two parabolae. This point has the same distance to two sites. Hence, it will lie on a Voronoi-edge. If it lies on three parabolae at the same time, then it has the same distance to three points. Hence, it must be a Voronoi-vertex. As the sweep-line moves down, the breakpoints will trace the edges of the Voronoi-diagram. Whenever three parabolae meet, the breakpoint becomes a Voronoi-vertex. This happens exactly, when the sweep line touches the bottom point of the circle through these three sites, because in this moment, any point which is to come further down will have a greater distance to the three sites than the three sites have among themselves. \ / x has the same distance to p and q. y has \ x q y / the same distance to q and r. Both trace a \ p /\___/\ r / Voronoi-edge. If x and y meet and q's arc \___/ \___/ disappears, a Voronoi-vertex is found. SL_______________________ Adjacent breakpoints: Two breakpoints b1 and b2, such there is no breakpoint with an x-coordinate between b1.x and b2.x. Arc of a beach line: The set of points on the beach line, which have an x-coordinate between the x-coordinates of two given adjacent breakpoints. Voronoi-sweeping: The following algorithm, which uses a beach line to calculate the Voronoi-diagram of a set of points: // This is only a very restricted algorithm (returns only vertices), // which has to involve plenty of pseudo-code due to its complexity. // Sorry for that... FUNCTION voronoi(P : Set of Point) : List of Vertex result: List of Vertex := <> // This structure stores the arcs of the beachline in the // leaves of a balanced binary tree. Such an arc is stored // in 'arc' by the corresponding site with x giving an // estimate of its x-position beachLine: Heap of (x:Real,arc:Point) with priority x := <> events: PriorityQueue of (p:Point,type:{SITE,VERTEX}, center:Point,arcIds:Set of Real) with the priority -p.y := <> // Insert all sites into the events-structure FOREACH p in P DO events.insert((p,SITE,nil,{})) // Run over first site (p,_,_) := event.deleteMin() beachline.insert((p.x,p)) WHILE !events.empty() DO (p,type,center) := events.deleteMin() IF type==SITE THEN beachLine.insert((p.x,p)) // Get my two neighbors pred:=beachLine.pred((p.x,p)) succ:=beachLine.succ((p.x,p)) // Find the closest arc among them IF |p-pred.arc|<|p-succ.arc| THEN // All vertex-events involving pred.arc must die... // Patch a piece of pred right to p beachLine.insert(( (p.x+succ.x)/2, pred) ELSE // All vertex-events involving succ.arc must die... // Patch a piece of succ left to p beachLine.insert(( (p.x+pred.x)/2, succ) // Find circles with predecessors and successors predpred:=beachLine.pred(pred) succsucc:=beachLine.succ(succ) disc:=smallestDisc3(predpred.arc,pred.arc,p) // Is this a good one? IF disc.center.y-disc.radius beachLine.delete(above) result := result \/ center RETURN result This algorithm works as follows (if it works at all in this version...): All sites are events. The beach line is kept in a heap (in the normal algorithm, it is kept in a balanced binary tree). When the beach line meets a site p2, it patches the beach line as follows: p1 p1 \____/ ~~~~~~> \_____/ \/ sweep line --------p2----- --------p2------- beachLine: ....p1.... ...p1,p2,p1.... Moreover, all vertex events which have been introduced on basis of the arc, which is just being split up, have to be removed: The arc can no longer contribute to any voronoi cell, because p2 is i its way. Furthermore, it is checked whether the sites of the neighboring arcs of the new arc make up a circle, which crosses the sweep line. If this is the case, then the bottom of this circle is added as a vertex event. If the sweep line reaches such a vertex event, then no further site below the sweep line can influence the area of the circle. Hence, its center is declared a vertex. Moreover, the arc directly above this vertex is deleted, as it contributes no longer to the beach line. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Delaunay-Triangulations ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Triangle: A polygon with 3 vertices. Quadrilateral: A polygon with 4 vertices. As for polygons, quadrilateral vertices are given in clockwise order. Skinny triangle: A triangle with one of its angles being small. Triangulation of a set of sites P: A maximal straight planar graph with * the sites as its vertices * line segments between the sites as edges * each bounded face being a triangle Every triangulation of P has the same number of triangles. Good triangulation: A triangulation with very few skinny triangles. Good triangulations improve interpolated function values, if only the function values of the vertices are known. Example: One knows the height of certain landscape locations and one would like to estimate the height of other locations. One can find a triangulation among the known locations and calculate the heights on the edges as weighted mean-values of their end points. This aproximation will be more realistic with good triangulations. Angle vector of a triangulation T: The tuple of all angles of all triangles of T, sorted in ascending order. One defines (a1,a2,...am) < (b1,b2,...bm) :<=> a1 < b1 || a1==b1 && (a2,...am)<(b2,...bm) Good triangulations have a large angle vector. Delaunay-Triangulation of a set of sites P: A triangulation of P with a maximum angle vector. The general position assumption will be made. The following properties hold: * The triangulation is a kind of inverse of the Voronoi-Diagram: The edges of the Delaunay-Triangulation connect exactly those points of P separated by one Voronoi-edge. * Three sites p1, p2, p3 are the vertices of a Delaunay-triangle iff the circle through them contains no site of P in its interior. * Two sites form a Delaunay-edge, iff there is a circle with the two sites on its boundary and no site in its interior. Notation: Del(P) Additional information: The convex hull can also be calculated for a set of 3-dimensional sites. Then the convex hull is a set of two-dimensional polygons, which fit together like the faces of a football. By the help of a 3-dimensional convex-hull-algorithm, one can calculate the Delaunay-Triangulation of a 2-dimensional set of points: Just imagine a paraboloid structure flying above the plane. This is as if someone placed a salad-bowl on the plane. Now have all 2-dimensional sites fly vertically upwards until they hit the salad-bowl. Calculate the convex hull of these 3-dimensional points. This gives a set of triangles between the sites, sticking on the salad-bowl. Now have these triangles fall back on the plane: This is the Delaunay-Triangulation. Proof: Consider the three sites of a triangle on the salad-bowl and draw a circle through them. There cannot be another site inside it, because all points interior to the circle have a greater distance to the salad-bowl center than the points on its boundary. Thus, it does not lie to the interior side of the triangle, which contradicts the principle of the convex hull. Illegal edge of a quadrilateral : A line segment p[i]--p[i+2], such the circle defined by p[i] and p[i+2] and a third vertex includes the fourth vertex in its interior. Flipping an edge in a quadrilateral : Deleting an edge p[i]--p[i+2] and inserting the edge p[i+1]--p[i+3]. /|\ /\ / | \ ~~~~~~~~> /__\ \ | / \ / \|/ \/ Flipping lemma: The fact that the Delaunay-triangulation is the triangulation, which does not contain any illegal edge. Flipper-algorithm: The following randomized incremental algorithm for the calculation of the Delaunay-Triangulation: FUNCTION checkEdge(a--b:LineSegment,p:Point, VAR triangles:Set of Triangle) // Does the edge a--b concern any other triangle? IF Ex c: in triangles THEN IF illegal(a--b,) THEN // Flip the illegal edge result := result \ {, } \/ {, } // Check the new edges checkEdge(a--c,p,result) checkEdge(b--c,p,result) FUNCTION flipper(P : Array of Point) : Set of Triangle // Add a dummy triangle, which bounds P p1 : Point := (min({p.x | p in P})-100, min({p.y | p in P})-100) p2 : Point := (p1.x (max({p.y | p in P})-p1.y)*2+p1.y+200) p3 : Point := ((max({p.x | p in P})-p1.x)*2+p1.x+200, p1.y) result : Set of Triangle := {} shuffle(P) FOR i:=0 TO P.length-1 DO // Determine the triangle which p[i] lies in := findTriangle(result, p[i]) // findTriangle can be constructed by // planar-point-location, see below // Split the triangle into 3 triangles around p[i] result := result \ {} \/ {, , } // Check and flip the new edges checkEdge(a--b,p[i],result) checkEdge(b--c,p[i],result) checkEdge(c--a,p[i],result) // Return result without dummy triangles RETURN result \ { | a in {p1,p2,p3} or b in {p1,p2,p3} or c in {p1,p2,p3} } This algorithm starts off with a large triangle, which bounds all sites. Then it adds one site after the other in random order: It determines the triangle a site lies in and splits this triangle into 3 new triangles around the site. Thereby, the edges of the old triangle might become illegal. We check each of them. Checking an edge a--b with respect to the new point p goes as follows: Find the triangle and find out whether the edge a--b is illegal in the quadrilateral . If this is the case, flip the edge a--b. Now, check each of the edges a--c and b--c, again with respect to p. 1) a________ is the triangle (sorry) which the /|\ \ new point p lies in. We added a--p and c/ | \p \d b--p. Now a--b is illegal in . \ | / / \|/_______/ b 2) a_______ We flipped a--b, but now a--c and b--c / \ \ might be illegal in their respective c/___\p \d quadrilaterals. Note that if a--c or b--c \ / / is flipped, the new edge will involve p. \ /______/ b The run-time of this algorithm is O(n*log(n)), if n is the number of sites. Locating the triangle which a site is in is possible in time log(n). This is done once for each site, resulting in time O(n*log(n)). When a new site is added, one edge may flip. After this flip, new flips can be necessary. Each of these flips adds an edge and each of these edges is incident to p. Thus, the degree of p tells us approximately how many flips have been done. The sum over the degrees of all nodes is a measure of the total number of flips. (I am getting hungry when I talk so much about Flips). The total number of edges is in O(n), such that all flipping also takes time O(n). ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Planar Point Location ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Planar subdivision: A straight planar graph, in which each node has at least degree 1. (?) Planar-Point-Location-Problem: The problem of, given a planar subdivision P and a point q, determining the face of P which q lies in. This is mostly done by preprocessing the planar subdivision in a convenient manner, so that afterwards, point location queries can be done fast. Slab between the real numbers x1 and x2: The set {(x,y) | x in [x1;x2]}. That is: The area between two vertical lines. Slab-Algorithm: An algorithm for solving the planar point location problem, which partitions the plane into slabs. Each point of the subdivision graph lies on the boundary of a slab left of it and on the boundary of a slab right to it. If a point q is to be located, its slab is found by binary search. Then, its face is determined by binary search within the slab. This algorithm can answer queries in time O(log(n)), but requires memory space in O(n^2), if n is the number of graph nodes. Parallel lines: Two lines with the same direction. Trapezoid: A quadrilateral, two boundary line segments of which are parallel. As quadrilaterals, trapezoids may have one boundary line segment of size 0. Shooting line of a vertex p in a planar subdivision P=(V,E): The set SL(p)={(p.x,y) | y in [p.y-a;p.y+b]}, with maximal a and b and so that |SL(p) /\ E'|==1, where E' is the set of all edge points in E. That is: A vertical line segment through p, which reaches up from p to the next edge of P and down from p to the next edge of P. Shooting trapezoid in a planar subdivision P: A trapezoid, which is given by the edges of P and the shooting lines of all vertices of P. Each shooting trapezoid is defined by 4 entities: * a line segment on the top * a line segment on the bottom * the vertex of the left shooting line * the vertex of the right shooting line For the TrapMap-algorithm, a shooting trapezoid will be represented as CLASS ShootingTrapezoid top : LineSegment bottom : LineSegment leftp : Point rightp : Point urneighbor : ShootingTrapezoid // upper right neighbor lrneighbor : ShootingTrapezoid // lower right neighbor leaf : TrapSearch // Pointer to where I'm stored | T7 | T6 | | | | | | | Here, a,b,c,d are vertices of P and a|____________|b | a--b, b--c, c--d, d--a are edges of P. |\ T1| |\ | All vertical lines are shooting lines. | \ | T2 | \ | The resulting trapezoids are T1,...T7 | \ | |T3\ | | d\|_______|___\|c | | | | | | | T4 | T5 | Neighbor of a shooting trapezoid A: Another shooting trapezoid, which has points of its boundary in common with A. Trapezoid search structure: A directed acyclic graph with one node called root, from which all nodes are reachable. The leaves of this graph (i.e. nodes of out-degree 0) contain shooting trapezoids. The other nodes are all of out-degree 2 and store either a shooting point ("x-node") or a line segment ("y-node"). CLASS TrapSearch type : {leaf, xnode, ynode} trap : ShootingTrapezoid // In case of leaf point : Point // In case of xnode line : LineSegment // In case of ynode left : TrapSearch right : TrapSearch FUNCTION locate(q:Point) : ShootingTrapezoid SWITCH this.type CASE xnode: IF q.x>this.point.x THEN RETURN this.left.locate(q) ELSE RETURN this.right.locate(q) CASE ynode: IF q.y above this.line THEN RETURN this.left.locate(q) ELSE RETURN this.right.locate(q) CASE leaf: RETURN this.trap This structure represents a trapezoid map. A point can be located by determining whether it is right or left of a shooting line or above or below a line segment. Depending on this decision, the search is continued in a sub-region, where there are again shooting lines and line segments. TrapMap-Algorithm: The following randomized incremental algorithm for the planar point location problem: FUNCTION makeTrapMap((V,E) : PlanarSubDivision) : TrapSearch // Add bounding box minx := min({p.x | p in V}) maxx := max({p.x | p in V}) miny := min({p.y | p in V}) maxy := max({p.y | p in V}) bb : ShootingTrapezoid := (top: y=maxy+100, bottom: y=miny-100, leftp : (minx-100,miny-100), rightp : (maxx+100,miny-100)) D : TrapSearch := (type: leaf, trap:bb) bb.leaf := D // Shuffle edges S : array of LineSegment := shuffle(array(E)) // Add edges one after the other FOR i:=0 TO S.length DO // Sort line end points from left to right IF S[i].p1.x>S[i].p2.x THEN swap(S[i].p1,S[i].p2) // Determine the trapezoid of S[i].p1 curTrap : ShootingTrapezoid := D.locate(S[i].p1) // Now split the trapezoid | _________________| This is how the trapezoid | | is being split into | above | three trapezoids , | | and |newLeft p1~~~~~~~~~~~~S[i] | | | below | |________________|_____ above : ShootingTrapezoid := (leftp:S[i].p1, top:curtrap.top, bottom:S[i]) below : ShootingTrapezoid := (leftp:S[i].p1, top:S[i], bottom:curtrap.bottom) newLeft : ShootingTrapezoid := (leftp:curtrap.leftp, rightp:S[i].p1, top:curtrap.top, bottom:curtrap.bottom, urneighbor:above, lrneighbor:below) curTrap.leaf:=(type: xnode, left: (type:leaf, trap:newLeft) right: (type:trap, trap:curtrap)) newLeft.leaf := curTrap.leaf.left above.leaf:=(type:leaf, trap:above) below.leaf:=(type:leaf, trap:below) // Move curtrap to be the imaginary trapezoid around // and curtrap.leaf:=curtrap.leaf.right // Replace all trapezoids along the line WHILE curTrap.rightp left of S[i].p2 DO // All trapezoids intersected become y-nodes. // They divide into and curtrap.leaf.type:= ynode curtrap.leaf.line:= S[i] curtrap.leaf.left:= above.leaf curtrap.leaf.right:= below.leaf IF curTrap.rightp above S[i] THEN // We are hit by a shooting line from above. // Thus, we must finish the current and // create a new one right besides above.rightp:=curtrap.rightp above.urneighbor := curtrap.urneighbor above.lrneighbor := (leftp:curtrap.rightp bottom:S[i], top:curtrap.lrneighbor.top) above := above.lrneighbor above.leaf:=(type:leaf, trap:above) curTrap:=curTrap.lrneighbor ELSE // We are hit by a shooting line from below below.rightp:=curtrap.rightp below.lrneighbor := curtrap.lrneighbor below.urneighbor := (leftp:curtrap.rightp bottom:curtrap.lrneighbor.bottom, top:S[i]) below:=below.urneighbor below.leaf:=(type:leaf, trap:below) curTrap:=curTrap.urneighbor // Now curTrap is the trapzoid which S[i].p2 lies in. // Split it into and and newRight : ShootingTrapezoid := (leftp:S[i].p1, rightp:curtrap.rightp, top:curtrap.top, bottom:curtrap.bottom, urneighbor:curtrap.urneighbor, lrneighbor:curtrap.ulneighbor) curTrap.leaf:=(type: xnode, left: (type:ynode, line:S[i], left:above.leaf right:below.leaf) right: (type:trap, trap:newRight)) RETURN D This algorithm maintains a Trapezoid Search Structure D. In the beginning, the structure just contains a large dummy trapezoid, which bounds the whole planar subdivision. Now, the edges are added incrementally. An edge intersects a number of trapezoids. D is already powerful enough to locate the trapezoid which contains the left starting point of the edge. This trapezoid is split into three trapezoids: One to the left of the starting point, one above the line and one below the line. Now, we follow the line. We replace each trapezoid we intersect by a y-node, which points to the new above trapezoid and the new below trapezoid. Whenever we are hit by a shooting line from above, we need to close the trapezoid above. When we are hit by a shooting line from below, we need to close the trapezoid below. When we reached the trapezoid which contains the line end point, we split this one similarly to the very first one: It becomes an x-node, the left side of which points to the above trapzeoid and the below trapezoid and the right side of which points to the remaining part of the original trapezoid right of the line end point. The expected query time is O(log(n)). Proof: We look at a path through D to locate a point q. Let the random variable X[i] denote the number of nodes created on this path by the i-th iteration of the algorithm. Then the expected path length is pathlength = E(SUM {X[i] | i=1..n}) = SUM {E(X[i]) | i=1..n} where n is the number of line segments. Any iteration i adds at most 2 nodes to the search path for q. This is the case when q lies in the trapezoid above or below the starting point of line S[i]: We first ask whether q is left or right of the starting point, and if it is right of the starting point, we need to ask whether q is above or below S[i]. These two additional case distinctions add 2 nodes. // Why does the script talk of 3 nodes (?) Anyway, it's just a // constant factor... Let P[i] be the probability that the i-th iteration contributes *any* node to the search path. Then we have E(X[i]) =< 2*P[i] The i-th iteration contributes a node to the search path of q exactly if a new trapezoid is added around q. This probability can be determined by backward analysis: Look at the trapezoidal map after the i-th iteration. We know that the line S[i] was added, but S[i] is any of the present lines with equal probability. Now take S[i] away to simulate what was before the i-th iteration. Did the trapezoid around q disappear? The probability of this trapezoid to disappear is equal to the probability that S[i] contributed a top, a bottom, a leftp or a rightp to this trapezoid. There are 4 lines involved in these data and there is an overall number of i lines, thus the probability of S[i] to contribute to the trapzoid is 4/i. This is the probability P[i]. Now we can calculate pathlength = SUM { E(X[i]) | i=1..n} =< SUM { 2*P[i] | i=1..n} =< SUM { 2*4/i | i=1..n} in O(log(n)) The expected size of D is in O(n). Proof: The number of leaves of D is the number of trapezoids. The number of trapezoids is in O(n), since every line adds only a constant number of trapezoids. Let k[i] be the number of trapezoids created in the i-th iteration. Then the number of inner nodes created in this iteration is k[i]-1. Since SUM {k[i] |i=1..n} in O(n), it follows that SUM {k[i]-1 |i=1..n} in O(n) and hence the overall number of nodes is in O(n). // Is it that easy (?) Euclidean Graph: A complete weighted undirected graph, the nodes of which are points and the edge weights of which are the euclidean distances of the nodes. EMST, Euclidean minimal spanning tree: A minimal spanning tree in an euclidean graph. The EMST is a subgraph of the Delaunay Triangulation of the vertices. Proof: We show that any edge (p,q) of the EMST is also an edge of the Delaunay Triangulation. We look at the smallest circle passing through p and q, i.e. the circle with p--q as its diameter. Now suppose there is a vertex r inside this circle. Then we have |(p,r)| < |(p,q)| and |(r,q)| < |(p,q)|. Due to the cycle property of MSTs, we know that the heaviest edge on a cycle is not needed for an MST. Hence, (p,q) does not belong to the MST, because it is the heaviest edge on the cycle p->q->r->p. This is a contradiction. Hence there cannot be a vertex inside the circle and (p,q) is a Delaunay edge. This yields an O(n*log(n))-algorithm for determining the EMST of n points: FUNCTION emst(P : Set of Point) : Set of Edge triangles := flipper(P); edges := {(u,v) | (u,v) is an edge in an element of triangles} RETURN kruskal((P,Edges,|.|)) We first calculate the Delaunay Triangulation and then the MST in it. Closest-Pairs-Problem: The problem of, given a set of sites, determining for each site its closest site. In the euclidean graph, this translates to finding for each vertex its shortest edge. These edges are part of the MST of this graph: Consider a vertex p[i] and its closest vertex p[j]. We can draw a cut around p[i] and the Minimum-Cut- Property tells us that the smallest edge leaving the cut must be part of the MST. This is p[i]--p[j]. Hence we can construct the EMST of the sites and traverse its edges to find the closest pairs: FUNCTION closestPairs(P : Array of Point) : Array of Point closestPoint : Array[P.length] of Point distance : Array[P.length] of Real := <+oo, +oo, ...> edges := emst(set(P)) FOREACH (p,q) in edges DO IF |(p,q)|1 DO step := step / 2 IF p right of P[i] THEN i := i + step ELSE i := i - step IF P[i+1] right of P[i] && p above P[i]--P[i+1] || P[i-1] left of P[i] && p above P[i-1]--P[i] THEN RETURN false i := rightmost step := rightmost-leftmost WHILE step>1 DO step := step / 2 IF p right of P[i] THEN i := i - step ELSE i := i + step RETURN ~(P[i-1] right of P[i] && p below P[i]--P[i+1] || P[i+1] left of P[i] && p below P[i-1]--P[i]) This algorithm does a binary search on the upper part of the polygon to find the polygon vertex with the closest x-coordinate. If the point lies above the corresponding edges, then it cannot lie inside the polygon. The same procedure is used for the lower part of the polygon. Largest-Empty-Circle-Problem: The problem of, given a set of sites, determining the largest circle with its center in the convex hull, such that no site lies in its interior. The center of this circle is either voronoi vertex or an intersction point of a voronoi edge and the convex hull. Proof: Suppose no site lies on the boundary of the circle. Then we can keep the center and enlarge it until one site is on it. Thus, this circle could not be maximal. Similarly, one can always enlarge circles with only 1 site on it or 2 sites on it. Thus, there must be 3 sites on the circle, which means that it is a voronoi vertex. Or there must be two sites on it and we are not allowed to enlarge the circle any more. The two sites on it mean that the circle center is on a voronoi edge. The fact that we are not allowed to enlarge it any more means that the center is about to cross the convex hull. Hence the center is on the intersection of a voronoi edge and the convex hull. The circle can be found as follows: FUNCTION LargestCircle(P : set of Point) : Circle ch : Polygon := convexHull(P) (V,E) : Graph := voronoi(P) // We need the whole graph here candidates : Set of Circle := {} // Make all interior voronoi vertices candidates FOREACH p in V DO IF p inside ch THEN candidates := candidates \/ {(p,)} // Find intersections of voronoi edges with convex hull intersections : Set of Point := lineIntersect(E \/ {ch[i]--ch[i+1]|i=0..ch.length}) FOREACH p in intersections DO candidates := candidates \/ {(p,)} result : Circle := ((0,0),0) FOREACH c in candidates DO IF c.radius > result.radius THEN result := c RETURN result All of the procedures called run in time O(n*log(n)), if n is the number of sites. All loops run in time O(n). Hence the algorithm runs in time O(n*log(n)). Post-office-problem: The problem of, given the location of n post-offices as points in the plane, determining the location for a new post-office, which has the greatest distance to all other post-offices, while at the same time lying inside the convex hull of the other post-offices. This problem corresponds to the Largest-Empty-Circle-Problem. The center of the largest empty circle fulfills the conditions for the new post-office. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Optimization ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Triangle inequality: The phenomenon that in a weighted graph (V,E,c), c(u,w) =< c(u,v) + c(v,w) for all u,v,w in V. Optimization problem: A tuple of a set L of "feasible solutions", a "cost function" f:L->Real and an element of {min,max}. The goal will be to find the feasible solution with the cost either maximized or minimized. Often, the set L of feasible solutions is not known explicitly, but only by some constraints which the elements of L fulfill. Many of the problems treated above can be seen as optimization problems: * Shortest paths: L is the set of possible paths, f is their length, to be minimized * Minimum spanning trees: L is the set of spanning trees, f is the sum of their edge weights, to be minimized * Maximum flow: L is the set of feasible flows, f is the flow value, to be maximized * Maximum cardinality matching: What's this (?) * Maximum weight matching: L is the set of all matchings, f is their weight, to be maximized * Linear programs: L is the set of feasible solutions and f(x)=x*c is to be maximized. Global optimum, Optimal solution of an optimization problem (L,f,MINMAX): An element x* of L, such that f(x*) = MINMAX({f(x) | x in L}). Minimization problem: An optimization problem, where the costs are to be minimized. Maximization problem: An optimization problem, where the costs are to be maximized. Every maximization problem can be translated to a minimization problem by negating the cost function -- and vice versa. Instance of an optimization problem: A certain optimization problem, with a given set of feasible solutions and a given cost function. Decision problem of an optimization problem (L,f,minmax): The problem of determining whether there is a feasible solution, which has a better cost than a given value ^x. That is Ex x in L: f(x) = minmax(f(x),^x) ? Oracle of an optimization problem (L,f,minmax): A function, which, given a real value b, returns * a feasible solution x with f(x) better than or equal to b, if there is any such solution * nil else. The oracle can be used to find the optimal solution, if f maps only to integers: FUNCTION optimalSolution((L,f,min):OptimizationProblem, P : FUNCTION Integer -> L) : L upperbound := 1 WHILE P(upperbound)==false DO upperbound := upperbound * 2 step : Integer := upperbound/2 WHILE step>0 DO IF P(upperbound)!=nil THEN upperbound := upperbound - step ELSE upperbound := upperbound + step step:=step/2 RETURN P(upperbound) This algorithm is designed for minimization problems. It first seeks an upper bound for the cost of the optimal solution. This is done by checking P for the powers of 2. As soon as the oracle P says that there is a feasible solution below any bound 2^k, we know that 2^k is an upper bound for the cost of the optimal solution. Then we determine the smallest x for which P still knows a smaller solution by binary search. This algorithm takes n*log(n) as much time as P, if n is the cost of the optimal solution. Integer Linear Program, ILP: A linear program with the additional constraint that the vector x be composed of integers. ILPs are NP-hard. Notational variant: x in {0,1,..k}^n := x[i] >= 0 for i=1..n x[i] =< k for i=1..n x[i] integer for i=1..n Multicommodity flow problem: The following maximization problem: Given * a weighted graph (V,E,cap) * k pairs of source node and target node (s[i],t[i]) in V x V find k flows f[i]:E -> Real+, such that * All e in E: All i in {1,...k}: f[i](e) >=0 * All e in E: SUM f[i](e) =< cap(e) * All i in {1,...k}: All v in V\{s[i],t[i]}: SUM {f[i](u,v) | (u,v) in E} = SUM {f[i](v,) | (v,w) in E} Maximize the overall outflow of the source nodes SUM {f[i](s[i],v) - f[i](w,s[i]) | i in {1,...k}, (s[i],v) in E, (w,s[i]) in E} The above formalization shows that this problem can be expressed as an LP. Set covering problem: The following minimization problem: Given * a set M={0,...m} * n subsets M[i] < M, such that M[1] \/ M[2] \/ ... \/ M[n] = M * a cost for each subset: c[1],...c[n] in Real+ select * a number of subsets M[i], such that their union is still M Minimize the number of required subsets. This problem can be formalized as an ILP: Minimize (c[1],...c[n])*(x[1],...x[n]) x[1],....x[n] in {0,1} define isin[i][j] := (i in M[j])?1:0 (isin[1][1],isin[1][2],...isin[1][n]) * (x[1],...x[n]) >= 1 ... (isin[m][1],isin[m][2],...isin[m][n]) * (x[1],...x[n]) >= 1 Here, x[i] is 1 if M[i] is selected and 0 else. Thus, the minimization term reflects exactly the sum of the costs of the selected subsets. Then we postulate that the x[i] are either 0 or 1. In order to ensure that each element of M appears in at least one of the selected subsets, we calculate the number of selected subsets in which the element occurs. This number should be greater than 1. Tutor-problem: The following maximization problem: Given * n students 1,2,...n * m tutors 1,2,...m * a happiness of each student for each tutor happy:{1,...n} x {1,...m} -> Real+ * a maximum number of students per tutor K Find * a mapping of students to tutors f:{1,...n}->{1,...m}, such that each tutor has at most K students All t=1...m: |{f(x)=t|x=1...n}| =< K Maximize the overall happiness, i.e. SUM {happy(s,f(s))|s=1...n}. This problem can be expressed as an ILP: We use n*m variables T[i][j] := 1 if student i is assigned to tutor j, 0 else We express the happy-function as an n*m-matrix H H[i][j]=happy(i,j) Then the problem is maximize SUM {H[i][j]*T[i][j] | i=1..n, j=1..m} T[i][j] integer for i=1...n, j=1..m SUM {T[i][j] | j=1..m} = 1 for i=1..n SUM {T[i][j] | i=1..n} =< K for j=1..m The sum in the first line is the sum of the total happiness under the current assignment T. The other two sums express that each student shall be assigned to exactly one tutor and that no tutor shall have more than K students. The problem can also be expressed as a min-cost-flow problem: We build a cost-capacity-graph with * the nodes V={s[1],...s[n]} \/ {t[1],...t[m]} \/ {target}, i.e. one for each student, one for each tutor and one target node * the edges E={s[1],...s[n]} x {t[1],...t[m]} \/ {t[1],...t[m]} x {target}, i.e. one from each student to each tutor and one from each tutor to the target node * the capacity cap:E->Real+ with * cap(s[i],t[j])=1 for i=1..n, j=1..m * cap(t[j],target)=K for j=1..m * the costs cost:V->Real with * cost(s[i],t[j]) = -happy(i,j) for i=1..n, j=1..m * cost(t[j],target)=0 * the supply supply:V->Real with supply(s[i])=1 for i=1..n supply(t[j])=0 for j=1..m supply(target) = -n Each feasible integer flow in this graph represents a valid mapping from students to tutors. The minimal integer cost flow is the one which maximizes the happiness. Traveling salesman problem, Traveling saleswoman problem, traveling salesperson problem, TSP: The following optimization problem: Given a weighted graph, find the shortest hamiltonian cycle in it. The idea is the following: Given a set of cities (the nodes) and connections between them (the links) with a certain distance (the weights), find the shortest way to visit all cities. Graph Coloring problem: The following optimization problem: Given a graph, find a mapping of nodes to colors, such that no two adjacent nodes get the same color. The goal is to minimize the number of required colors. Mixed Integer Linear Program, MILP: A linear program with the additional constraint that some of the components of the vector x be integers. MILPs are NP-hard. Linear relaxation of a ILP: The ILP without the constraints that the vector x be composed of integers. Knapsack-problem: The following optimization problem: Given a set of objects, each with a real weight and a real profit, determine a subset of them, such that the sum of their weights does not exceed a given limit. The goal is to maximize the sum of the profits. The knapsack- problem can be seen as an ILP: Maximize (profit[1], profit[2],...profit[n])*x (weight[1], weight[2],...weight[n])*x =< maxweight x in {0,1}^n The components of x indicate whether the i-th object is to be selected (x[i]=1) or not (x[i]=0). __________________ / ? V _ _ | | | | _ | | | | | | | | | | _ | | <- knapsack |_| |_| |_| |_| |___| 10$ 8$ 12$ 2$ Profit-density of a knapsack object: Its profit divided by its weight. Fractional knapsack-problem: The linear relaxation of the knapsack-problem. That is: Objects may also be taken in fractions. This problem can be solved by the following algorithm: FUNCTION linfracKnapsack(object: Set of n (profit,weight), maxweight : Integer) : Set of (profit,weight,fraction:Real) sorted := sortBy(object,profit/weight,DECREASING) result : Set of (profit,weight,fraction:Real) w := 0 FOR i:=0 TO n-1 DO IF sorted[i].weight+w < maxweight THEN result := result \/ {(sorted[i],1)} w := w + result[i].weight ELSE result := result \/ {(sorted[i],sorted[i].weight/w)} BREAK RETURN result This algorithm first sorts the objects according to their profit-density. Then it adds them one after the other until the maximum weight is bypassed. Then the last object is added only partially. The running time is in O(n*log(n)). If the partial object is left away, the resulting knapsack composition is not necessarily an optimal solution of the original problem. But the profit achieved in linear relaxation is an upper bound for the real knapsack problem. D ________________ C _ / V B _ | | | 1 A | A _ | | | | | 1 B | _ | | | | | | | 2/3 C | <- knapsack |_| |_| |_| |_| |_______| 10$ 15$ 20$ 20$ 1kg 2kg 3kg 4kg 5kg The problem can also be solved in expected time O(n) by the following algorithm: FUNCTION linfracKnapsack(object: Set of n (profit,weight:Real), maxweight : Integer) : Set of (profit,weight,fraction:Real) // Pick a Pivot-element pivot : (profit:Real,weight:Real) := eps(object) // Distribute into 3 sets greater : Set of (profit:Real,weight:Real) := {} greaterSum : Real := 0 equal : Set of (profit:Real,weight:Real) := {} equalSum : Real := 0 smaller : Set of (profit:Real,weight:Real) := {} FOREACH o in object DO IF o.profit/o.weight>pivot.profit/pivot.weight THEN greaterSum := greaterSum + o.weight greater := greater \/ {o} IF o.profit/o.weight==pivot.profit/pivot.weight THEN equalSum := equalSum + o.weight equal := equal \/ {o} IF o.profit/o.weight contains too many elements, recurse IF greaterSum > maxweight THEN RETURN linFracKnapsack(greater,maxweight) // Can we fill the knapsack with the elements from ? maxweight := maxweight - greaterSum WHILE equal != {} DO x := pickAny(equal) IF maxweight=< x.weight THEN RETURN (greater x {1}) \/ {(x,x.weight/maxweight)} greater := greater \/ {x} maxweight := maxweight - x.weight // Come here if we need to consider the elements from RETURN (greater x {1}) \/ linKnapsack(smaller,maxweight) This algorithm picks a pivot-element at random among the objects. Then it runs through the set of objects and partitions it into 3 sets S1, S2 and S3, according to whether the profit-density is better than, equal to or worse than the profit-density of the pivot-element. If the weight of the better set S1 exceeds the maximum weight, then the algorithm recurses on this set. Else, it tries to enlarge the set by adding objects from S2. If the knapsack is filled, it is returned. Else, the algorithm recurses on the remaining objects with a reduced maximum weight. This algorithm works because we get an optimal fractional knapsack if we chose the objects with large profit-densities first. It has an expected runtime of O(n). Proof: If the pivot-element is chosen truly at random, the sets S1 and S3 will have approximately the same size. Hence, the problem size is divided by two in each recursion. There may be at most log(n) recursions, where the i-th recursion takes O(2^-i*n) steps. Thus, the overall time is SUM {O(2^-i*n) | i=0...log(n)} in O(n). Brute-Force: The following blueprint for an algorithm for an optimization problem: FUNCTION bruteForce((L,f) : MinimizationProblem) : L // Start with a start state cand : Set of Set of L := INVARIANT Ex M in cand: x* in M, x* is solution of (L,f) bestsofar : L := nil WHILE cand != {} DO // Deleting curr : Set of L := cand := cand \ {curr} // Pruning (Bounding) curr := curr \ // Selecting IF |curr| == 1 THEN IF f(eps(curr)) better than f(bestsofar) THEN bestsofar := eps(curr) IF |curr| > 1 THEN // Splitting curr := curr[1] \/ ... \/ curr[k] such that curr[i]!=curr // Branching cand := cand \/ {curr[1],...curr[k]} RETURN bestsofar This algorithm maintains a set of sets of feasible solutions. It starts with the singleton set of any feasible solution. Then it always picks a set out of . Out of this set , obviouly useless feasible solutions are eliminated. If ony one feasible solution remains, this solution is compared to the best solution, which has been found so far. If consists of more than one feasible solution, then is split up into subsets. Each of these subsets is added as a new element to . This continues until no more candidates are left. Then, the best solution found so far is returned. This blueprint leaves much freedom: * the representation of L can be chosen * the start element of can be chosen * it can be chosen which elements of are deleted * it can be chosen which elements of are considered useless * it can be chosen how to split // Problems with discrete feasible solutions can also be treated as // graphs: The nodes are the feasible solutions. The links represent // the creation of a new feasible solution out of another one. (The // nodes can be calculated at runtime.) The source node is one known // feasible solution. All nodes with a good value for f are accepted // as target nodes. Now, a path from the source node to a target // node needs to be found. To do so, one can for instance use A*. // For a detailed treatment of this approach, see "Intelligence // artificielle" (ia.txt, in French) or "Artificial intelligence" // (ai.txt) Iterative deepening: The generic brute force algorithm, which is run with an additional parameter d:Integer. With each element of cand, one stores a integer value indicating how often this element has been subject to a split operation. Pruning is done as follows: If an element of cand has already been split d times or if it is known that the element will have to be split all in all more than d times, then this element is left away. Iterative deepening is first run with d=1. If it does not find a solution, it is run with d=2 etc. This corresponds to a breadth-first-search, which is simulated by repeated depth-first-searches. 15-puzzle: The following optimization problem: Given a 4x4 matrix with the numbers 0..16 in it, determine a sequence of swappings of the 0 with a neighboring element, such that the matrix takes the form 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 The goal is to find a minimum number of such swappings. This problem can be solved by iterative deepening as follows: * A feasible solution is a sequence of swappings, together with the resulting matrix. Hence we know for each feasible solution, how many swappings have already been made. * f is the number of elements not in the right place * is a stack of sets of feasible solutions * starts with the singleton set of the singleton set containing the feasible solution with 0 swappings * is a set of feasible solutions Since is a stack, always contains the newly expanded feasible solutions * Pruning is done by discarding feasible solutions, which will need more swappings than the parameter d of the iterative deepening call * Splitting is done by adding to each feasible solution of one more swapping We run the iterative deepening until a feasible solution has been found where all elements are in the right place, i.e. where f returns 0. // For a detailed treatment of the 15-puzzle, see "Intelligence // artificielle" (ia.txt, in French) ILP-Branch-and-Cut: The brute-force-algorithm applied to ILPs. This can be done as follows: * A feasible solution is a constraint. It represents the set of vectors fulfilling the constraint. * f is the cost function of the ILP * is a stack of sets of feasible solutions * starts with the singleton set, which contains the set of the constraints of the ILP * is a set of constraints. It represents the set of those vectors which satisfy all of the constraints. Since is a stack, always contains the newly expanded feasible solutions. * is a feasible solution of the ILP. * Pruning and selecting is done as follows: Find an optimal solution x' for the linear relaxation of . If f(x') is worse than f(bestsofar), then delete completely. If x' is composed of integers, then store x' in . * Splitting is done as follows: Consider a non-integer component x'[i] of x'. Produce two new sets of feasible solutions: curr \/ {x[i]==x'[i]} // I have the impression that 15-puzzles and ILPs cannot be solved // coherently by the same brute-force template. Types mismatch, some // steps are omitted (e.g. IF |curr|==1 THEN ...) and the two // approaches seem rather different. (?) Example: Solving a knapsack-problem by branch-and-cut maximize (10,15,20,20)*x (1,2,3,4)*x =< 5 0 =< x[i] =< 1 x integral cand : Stack of Set of Constraint := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1}> bestsofar : Vector of Integer:= nil // cand !={} ~~~~> enter loop // pick first element of cand curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1} cand := <> // Find optimal solution for linear relaxation of x' : Vector := (1,1,2/3,0) // Select and prune // f(x')=38 is better than f(bestsofar)=+oo // x' is not integral, it cannot be stored in bestsofar // Splitting on x[2]=2/3 curr1 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] =< 2/3} curr2 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] >= 2/3} // This simplifies immediately to curr1 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0} curr2 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} // Add and to cand := < {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0}, {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} > // cand !={} ~~~~> enter loop // pick first element of cand curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2]==0} cand := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1}> // Find optimal solution for linear relaxation of x' := (1,1,0,1/2) // Select and prune // f(x')=35 is better than f(bestsofar)=+oo // x' is not integral, it cannot be stored in bestsofar // Splitting on x[3]=1/2 curr1 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=0} curr2 := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=1} // Add and to cand := < {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=0}, {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=1}, {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} > curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=0}, cand := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=1}, {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} > x' := (1,1,0,0) // f(x')=25 better than f(bestsofar)=+oo // x' is integral ~~~~> store x' as new best solution bestsofar := (1,1,0,0) // No splitting is done curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 0, x[3]=1} cand := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} > x' := (1,0,0,1) // f(x')=30 better than f(bestsofar)=25 bestsofar := (1,0,0,1) // No splitting is done curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[2] == 1} cand := <> x' := (1,0,0,1) // f(x')=37 better than f(bestsofar)=30, but not integral // Split on x[1]=1/2 cand := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[1]=1, x[2] == 1}, {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[1]=0, x[2] == 1}> curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[1]=1, x[2] == 1} cand := <{(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[1]=0, x[2] == 1}> x' := (0,1,1,0) // f(x')=35 better than f(bestsofar)=30 bestsofar := (0,1,1,0) // No splitting is done curr := {(1,2,3,4)*x =< 5, 0 =< x[i] =< 1, x[1]=0, x[2] == 1} cand := <> x' := (1,0,1,1/4) // f(x')=35 is equal to f(bestsofar)=35, so skip x' // cand is empty, bestsofar=(0,1,1,0) Composable optimization problem: An optimization problem, the optimal solution of which can be seen as the composition of the optimal solution of a sub-problem and some additional element. In this case, the problem is said to obey the "principle of optimality". Examples: * The shortest path from graph node u to graph node v via a graph node w can be seen as a composition of the shortest path from u to w and the shortest path from w to v. * The optimal solution for a knapsack problem of weight maximum maxw can be seen as a composition of a knapsack problem with maximum weight maxw-w(x) plus the addition of the object x, where x is an object contained in the optimal solution. However, not all problems have compositionable optima: In the graph coloring problem for instance, the addition of a node makes the solution of the previous graph useless. Dynamic programming: The following blueprint for an algorithm for a composable optimization problem: Prerequisites: * Find integer measures i,j,k,... for the size of a sub-problem instance, such that the instance is trivially solvable for i=k=...=0 and such that i=upper1,j=upper2,... corresponds to the original problem. Example: For the knapsack-problem, measures may be * the number of objects available * the maximum weight of the knapsack If all objects are available and the weight is equal to the maximum weight of the original knapsack, then the sub-problem corresponds to the original problem. * You may find additional parameters, which do not influence the problem size, but characterize the solution of a sub-problem in a way. In the following, the term "measures" will refer to these parameters as well. * Introduce a function w, which maps the integer measures to the optimal cost of the sub-problem instance described by the integer measures. This function should be trivial for i=j=...=0. For other values, the function should be defined recursively. The recursion may choose the smaller sub-problem, on the basis of which the current w is calculated. Example: For the knapsack-problem, it holds w(numobj,maxweight) = (numobj==0 || maxweight==0)?0: max(w(numobj-1,maxweight), w(numobj-1,maxweight-weight[numobj-1])+profit[numobj-1]) This means: The profit for the sub-problem with the first objects available and a maximum weight of can be calculated by the following choice: Either the -th object is left out of the knapsack. Then w(numobj,maxweight)=w(numobj-1,maxweight) Or the -th object is included in the knapsack. Then the profit is the optimal solution for the sub-problem with weight minus the weight of the -th object. The profit is then increased by the profit of the -th object. The algorithm is: FUNCTION dynamicProgramming((L,f,minmax):OptimizationProblem, m : array[n] of Integermeasure, upper: array[n] of Integer, w : FUNCTION Integer^n->Real) : L // Set of measures of all possible subproblems subproblems := {0..upper[0]} x {0..upper[1]} x ...{0..upper[n-1]} // Allocate an n-dimensional table profit : array[upper[0]]...[upper[n-1]] of Real := <0...> choice : array[upper[0]]...[upper[n-1]] of subproblems := // Fill profit table for trivial values FOREACH x in subproblems: w(x) can be computed trivially DO profit[x] := w(x) // Build up solutions to sub-problems incremetally FOR i[0]:=0 TO upper[0] DO FOR i[1]:=0 TO upper[1] DO ... FOR i[n-1]:=0 TO upper[n-1] DO // Calculate w(i[0],i[1],...i[n-1]) by using table entries profit[i[0]][i[1]]...[i[n-1]] := w(i[0],i[1],...i[n-1]) choice[i[0]][i[1]]...[i[n-1]] := position : Vector of Integer := (upper[0],upper[1],...upper[n-1]) result := nil WHILE w(position) is not trivial DO by regarding and choice[position]> position := choice[position] RETURN result This algorithm computes the optimal solution by starting with a trivial sub-problem instance. Then, the recursive definition of w is used to build up a table, which stores for each possible sub-problem the optimal cost. Together with the optimal cost, one stores the smaller sub-problem, on the basis of which this optimal cost has been calculated. When all sub-problems have been calculated, one looks at the sub-problem, which corresponds to the original problem. By help of the choice taken to calculate its cost, one can find on which smaller sub-problem the original sub-problem bases. This identifies a part of the solution. One goes on, tracing back in the table, until the whole solution has been constructed. Depending on the structure of the recursion of w, it can be possible to use a profit-table with less dimensions. It is also possible to write a recursive algorithm basing only on the recurrence of w. However, this would cause the function to be called again and again during the process with the same parameter values. The dynamic programming approach avoids this: It stores the results in the table, so that multiple future recurrence equations can use them. // For a real-world-application, see the Viterbi parser in // "Computational Linguistics" (cl.txt), which determines the // word class (like noun, adjective etc.) for the words of a text // by dynamic programming. Picture for knapsack: maxweight 0 1 2 3 ... C numobj 0 0 0 0 0 ... 0 <- trivial solutions 1 * * * * ... * <- computed values 2 * * * * ... * 3 * a * * b * '---weight(4)--| <- c = max(a,b) V 4 * * * * c ... n Dynamic programming solution for the knapsack-problem: Dynamic programming with the measures and w-function described above for a knapsack problem with integer weights. The basic assumption is that in the non-trivial case, w(numobj,maxweight) = max(w(numobj-1,maxweight), w(numobj-1,maxweight-weight[numobj-1]) +profit[numobj-1]) Proof: Consider the knapsack-problem with the first objects given and with a maximum weight of . The profit cannot decrease, if an additional object is available, because we can always keep the old composition of the knapsack. Hence w(numobj,maxweight) >= w(numobj-1,maxweight) If w(numobj,maxweight) = w(numobj-1,maxweight), we're done. Now assume w(numobj,maxweight) > w(numobj-1,maxweight) This means that the profit increases by making available the -th object. Hence this object must be present in the knapsack. Now assume that w(numobj,maxweight) < w(numobj-1,maxweight-weight[numobj-1]) +profit[numobj-1] Then we already know a better solution: Take the solution of w(numobj-1,maxweight-weight[numobj-1]) and add the -th object. The result will be w(numobj,maxweight) = w(numobj-1,maxweight-weight[numobj-1]) +profit[numobj-1] Now suppose that adding the -th object somehow resulted in a magic solution, which is better than what we predicted: w(numobj,maxweight) > w(numobj-1,maxweight-weight[numobj-1]) +profit[numobj-1] Then take the -th object out of this magic solution. The result is a profit of w(numobj,maxweight)-profit[-1] which, according to the assumption, is greater than w' = w(numobj-1,maxweight-weight[numobj-1]) But then w' cannot be optimal, because we just showed how to derive a better solution for -1 objects -- namely from the magic solution. Hence the magic solution cannot exist and the recursion equivalence of w results. The dynamic programming algorithm can be optimized as follows: FUNCTION knapsack(profit : array[n] of Real, weight : array[n] of Integer, maxweight : Integer) : Set of Integer // Sub-problems for a fixed number of objects // and for all possible weights table : array[maxweight] of Real := <0,0,...0> // choice[i][C]=present <=> in the sub-problem with i objects // and a maximum weight C, the i-th object is in the knapsack choice : array[n][maxweight] of {present,notPresent} := // Incrementally add objects FOR i:=0 TO n-1 DO // Now i objects are given. Compute the optimal weight for // each possible upper bound. FOR C:=maxweight DOWNTO weight[i] DO IF table[C-weight[i]] + profit[i] > table[C] THEN table[C] := table[C-weight[i]] + profit[i] choice[i][C] := present // Reconstruct the solution result : Set of Integer := {} C := maxweight FOR i := n-1 DOWNTO 0 DO IF choice[i][C]==present THEN result := result \/ {i} C := C - weight[i] RETURN result This algorithm makes use of the fact that the optimal profit of a sub-problem can be calculated basing on the sub-problem with one object less. Hence, only one line is necessary of the 2-dimensional array of the normal dynamic programming algorithm. The time is in O(n*maxweight), the space consumption is in O(maxweight+n) plus n*maxweight bits. MASS-Problem, Maximum Alternating Sum Subsequence-Problem: The following maximization problem: Given an input sequence of n numbers S=, find a subsequence s[i1],s[i2],... with i1 < i2 < ..., such that the alternating sum s[i1]-s[i2]+s[i3]... is maximized. This problem can be solved by dynamic programming. As a measure for the size of a sub-problem, we choose the length of the input sequence. We use two recursive functions, odd:{1,...n}->Real+ even:{1,...n}->Real+. odd(i) is the value of the maximum odd-length alternating sum subsequence of the first i elements of S. even(i) is the value of the maximum even-length alternating sum subsequence of the first i elements of S. We have odd(0) = -oo even(0) = 0 That is: If we have 0 elements, we cannot have a subsequence of odd length. With 0 elements, the value of the only even subsequence is 0. We can decribe the two functions recursively: odd(i) = max(odd(i-1),even(i-1)-s[i]) even(i) = max(even(i-1),odd(i-1)+s[i]) That is: In order to find the maximum odd subsequence in the first i elements, we can just leave s[i] out and thus have imply odd(i-1). Otherwise, we can choose to find the maximum even subsequence within the first i-1 elements and add the i-th element. Analogously, we can determine the value for even(i). These two functions can be reduced to one by an additional parameter: w(0,1) = -oo w(0,0) = 0 w(i,1) = max(w(i-1,1),w(i-1,0)-s[i]) w(i,0) = max(w(i-1,0),w(i-1,1)+s[i]) Now, we have a function w:Natural^2->Real+, which can be plugged into the dynamic programming algorithm. // Gernot Gebhard suggested that a greedy algorithm might also work. Greedy strategy: The following strategy for solving a problem: Interpret the solution of the problem as a sequence of decisions. Always take a decision that locally seems to be the best. Repeat this until the solution is complete, never undo a step. This strategy is employed very successful in the following algorithms: * Kruskal: The algorithm grows a spanning tree by always picking the best (i.e. shortest) edge * Jarnik-Prim: The algorithm also grows a spanning tree by expanding it always by the seemingly best edge * Dijkstra: The algorithm finds a shortest path by adding up the shortest edges * Fractional knapsack: This algorithm always chooses the object with the largest profit-density Often, the greedy strategy can be implemented by a polynomial time algorithm. However, it does not work for all problems and especially often not for the NP-hard ones. Example: Consider the knapsack-problem with object profit weight profit/weight Greedy Optimal 1 10 1 10 * 2 15 2 7.5 * * 3 20 3 6.7 * 4 20 4 5 ------------------ Profit: 25 35 Change-problem: The following minimization problem: Given * a set of coin values c[1],...c[n] >0 in Natural * an amount m in Natural Find * a set S of coin values, such that SUM S = m Minimize the number of coins in S, i.e. minimize |S|. For the usual coin values 1,2,5,10,20,50,100,200, a greedy algorithm works fine: FUNCTION greedyChange(c : array[n] of Natural,m : Natural) : MultiSet of c[i] c := sort(c,DECREASING) result : MultiSet of c[i] := {} FOR i:=0 TO n-1 DO WHILE m>c[i] DO result := result \/ c[i] m := m-c[i] RETURN result But if we allow a coin c[i]=4, then this does not work out: For an amount of m, the algorithm would return {5,2,1}, whereas the optimal solution would be {4,4}. The general problem can be solved by the following dynamic program: As a measure, we chose the amount m. The w-function is w(0) := 0 w(m) := min({w(m-c[i]) | i=0...n-1, m-c[i]>=0})+1 For an amount of 0, we need 0 coins. For an amount of m, we check for each coin i, whether m-c[i] could be assembled efficiently. We choose the coin i, for which m-c[i] needed the least number of coins. Then we add the coin i (i.e. the number 1). This recurrence can be plugged into the normal dynamic programming algorithm. CD-problem: The following maximization problem: Given * n songs 1,2,3,...n * a compatibility for each pair of songs compat:{1,...n}^2->Real+ This compatibility compat(x,y) measures how good the song y sounds if it is heard after song x. Find * a sequence of the n songs to put on a CD Maximize the compatibility within the sequence. This problem can be solved by dynamic programming: As a measure for the size of a sub-problem, we choose the subset of songs of this sub-problem. There are 2^n possible subsets of {1,...n}. (The subsets can e.g. be expressed as integers by writing a '1' for each selected song and a '0' for each unselected song and then reading the resulting sequence as a binary number.) As an additional parameter for a sub-problem, we choose the number of the first song on the CD. Then the recurrence is as follows: w({},x) = 0 w(s,x) = (x in s)? max({w(s\{x},z)+compat(x,z) | z=1...n}): -oo That is: If the subset of songs on the CD is empty, then the maximum compatibility achievable is 0. Now suppose we are given a non-empty subset s of songs and a song x, which shall be the first song in the CD containing s. What is the maximum compatibility we can achieve under these conditions? At first, x needs to be in s. Then, we check out all CDs with the songs s\{x}. We have every song in s\{x} be the first song on the CD and we look whether song x would sound good before this first song of s\{x}. This gives the maximal achievable compatibility for this sub-problem. The recurrence can be plugged into the dynamic programming algorithm. The final result has to be calculated afterwards as max({w({1,2,...n},x) | x=1...n}) The running time is O(n^2*2^n), because the table has a size of n*2^n and for calculating one entry, we need to scan all n songs. Euler in a Tree: An euler circuit in an undirected tree. It can be calculated by plugging the following DFS-Action into the DFS-Algorithm: CLASS eulerInTRee EXTENDS DFSAction ec : List of Edge := <> PROCEDURE traverse((u,v) : Edge) ec.pushBack((u,v)) PROCEDURE backtrack((u,v) : Edge) ec.pushBack((v,u)) Approximation-ratio of an algorithm for a minimization problem: The value r, such that for all problem instances I f(x(I)) / f(x*(I)) < r where x(I) is the solution produced by the algorithm and x*(I) is the optimal solution. r-approximation of a minimization-problem: An algorithm for the problem with an approximation ratio of r. 2-approximation for the TSP: The following algorithm for solving the TSP in an undirected graph (V,E) with triangle inequality with an approximation-ratio of 2: FUNCTION tspApprox2((V,E,c) : WUCGraph) : List of E mst : Set of E := kruskal((V,E,c)) ec := eulerInTree() dfs((V,E,c), eps(V), ec) // Cut out doubly visited nodes changed := true WHILE changed DO changed := false vistedNodes : Set of Node := {} ecelement : List of Edge := ec // This is a pointer to ec! WHILE ecelement!=<> DO vistedNodes := visitedNodes \/ {ecelement.car.startNode} IF ecelement.car.endNode in visitedNodes && ecelement.car.endNode == ecelement.cdr.car.startNode THEN changed := true ecelement.car.endNode := ecelement.cdr.car.endNode ecelement.cdr := ecelement.cdr.cdr RETURN ec This algorithm supposes non-negative edge weights. It first calculates an MST in the graph. Then, it considers the bidirected graph of this MST and finds an euler circuit in it. Out of this Euler circuit, it glues together repeatedly edges of the form (u,v) (v,w) to an edge (u,w), if v has been visited before. The circuit returned by this algorithm is a 2-approximation. Proof: Consider the optimal TSP-tour C'. If one takes out one edge, one gets a spanning tree T'. Since the edge weights are supposed to be positive, it holds that c(T')=v->w) =< c(u->w). If C* is the resulting tour, then c(C*)= RETRN EC This algorithm first finds a euclidean minimal spanning tree. Then, it collects all points with odd degree in this EMST. For these points, it finds a minimum euclidean weight matching. This is possible because the number of nodes with odd degree in a graph is even. Now the algorithm considers all edges of the EMST and all edges of the matching and finds an euler cycle in this graph. This is possible because all nodes have even degree: Every node with a formerly odd degree has one additional edge from the matching. Then, the algorithm shortcuts this cycle by skipping already visited nodes. This algorithm is a 1.5-approximation of the TSP. Proof: According to the proof used in the 2-approximation, the EMST is a sub-graph of the optimal TSP-tour C*, hence c(emst){0,...m-1} x : FUNCTION {0,...n-1}->{0,...m-1} := {} // Load of each machine L : Priority Queue of (load:Real,machine:Integer) with priority load := <> FOR i:=0 to m-1 DO L.insert((0,i)) // Jobs still to do J : Priority Queue of (time:Real,job:Integer) with priority -time := <> FOR i:=0 to n-1 DO J.insert((t[i],i)) WHILE J != <> DO nextjob := J.deleteMin() nextmachine := L.deleteMin() nextmachine.load := nextmachine.load + nextjob.time x(nextmachine.machine) := nextjob.job L.insert(nextmachine) RETURN x This is a greedy algorithm: It always picks the biggest non-assigned job and assigns it to the machine with the least load. The algorithm runs in time O(n*log(n)), but does not always yield an optimal solution. Makespan-lemma: The fact that the makespan of a scheduling problem can be bounded by Lmax =< SUM { t[j]/m | j=0..n-1 } + (m-1)/m*t[last] where is the index of the job finishing last. Proof: We know that if the job had been mapped to another machine, then this machine would take longer than x(last): All i!=x(last): L(i)+t[last] >= Lmax && L(x(last)) = Lmax Now, we sum up all L(i)+t[last] and L(x(last)) and get SUM { L(i)+t[last] | i!=x(last) } + L(x(last)) >= m*Lmax This is nothing else than the overall time of all jobs, plus (m-1) times t[last]: SUM { t[j] } + (m-1)*t[last] >= m*Lmax => SUM { t[j]/m } + (m-1)/m*t[last] >= Lmax |-----m------| _ _ |L| |L| _ |_| |_| |L| SUM{t[j]}+ (m-1)*t[last] >= m*Lmax _ _ |_| |_| |_| _ |_| |_| |_| |_| |_| |_| ^---- This is the machine x(last) ^----^--------- Other machines, with a supposed additional Two-job-LPT-lemma: The fact that if the optimal solution of a scheduling problem assigns at most two jobs to one machine, then LPT computes an optimal solution for this problem. Proof: Suppose the jobs are ordered by decreasing time. Then the LPT will take the first job and assign it to the first machine. It will take the second job and assign it to the second machine etc. This continues to the n/2-th job, which is assigned to the last machine. Then, the next job is also assigned to the last machine, the next job is assigned to the but-last machine etc. Jobs 1,2,3,4,5,6 _ |_|6 _ Machines A, B, C | | |_|5 _ | | | | |_|4 |_|1 |_|2 |_|3 A B C Now consider any optimal solution. This optimal solution can be transformed as follows whithout changing its optimality (suppose the variables to be the weights of objects, placed as the picture suggests): * a c ~~~~~> c a if a+b < c+d b ... d d ... b && d < b * a c ~~~~~> c a if a+b < c+d b ... d b ... d && d > b * c ~~~~~> c if d > b b ... d b ... d The result is the solution produced by the LPT. (Really (?)) LPT-approximation-lemma: The fact that LPT has an approximation-ratio of 4/3 - 1/3/m. Proof: Assume that the approximation ratio is worse for an LPT-solution x and a optimal solution x*, i.e. 4/3 - 1/3/m < f(x)/f(x*) = Lmax/f(x*) =< SUM {t[j]/m} /f(x*) + (m-1)/m*t[last]/f(x*) by Makspan-Lemma =< f(x*) /f(x*) + (m-1)/m*t[last]/f(x*) since f(x*) >= SUM{t[j]/m} =< 1 + (m-1)/m*t[last]/f(x*) 1/3 - 1/3/m =< (m-1)/m*t[last]/f(x*) 1/3*m - 1/3 =< (m-1)*t[last]/f(x*) 1/3 =< t[last]/f(x*) f(x*) =< 3*t[last] Since is the last job added by LPT, it is the shortest job. Since the machine which runs longest can only treat 3 times the shortest job, it follows that every machine has at most 2 jobs. But then, LPT is optimal according to the two-job-LPT-lemma. Hence f(x) = f(x*), which leads to the contradiction 4/3 - 1/3/m < f(x)/f(x*) 4/3 - 1/3/m < f(x*)/f(x*) 4/3 - 1/3/m < 1 1/3 - 1/3/m < 0 m < 1 r-Dominating set, r-pareto-optimal solution in a knapsack problem: The optimal profit achievable if the knapsack has a maximum weight of r. (?) Nemhauser-Ullmann algorithm: The following modification of the dynamic programming algorithm for the knapsack problem: Instead of calculating a two-dimensional table with the optimal profit for each possible number of available objects and each possible maximum weight, compute for each number of available objects a profit list. This profit list contains pairs of a maximum weight and the achievable profit: numobj 0 <(maxw:0,profit:0)> 1 <(maxw:0,profit:0), (maxw:w[0],profit:p[0])> 2 <(maxw:0,profit:0), (maxw:w[0],profit:p[0]), (maxw:w[1],profit:p[1]),(maxw:w[1]+w[2],profit:p[1]+p[2])> ... ... This is a kind of shortcut of the dynamic programming table. The list for the i-th object can be calculated as follows: Copy the preceding list up to the element with maxw>=w[i]. Then run through the remaining elements. Be L the current car/cdr-pair. Be L' the latest car/cdr-pair with L'.car.maxw+w[i]=0 and (f(x)-f(x'))/T>0. Hence exp((f(x)-f(x'))/T)>1 and x' is chosen with probability 1. Now suppose f(x')>f(x). For large T, (f(x)-f(x'))/T is close to 0, so that exp((f(x)-f(x'))/T) is large and x' has a good chance of being chosen -- although it is sub-optimal. But as T decreases, exp((f(x)-f(x'))/T) will become small and it is less likely that x' is chosen. As a result, the algorithm is very curious in the beginning. It explores new feasible solutions, although they are bad. In the end, the algorithm becomes more eager for optimization and is no longer willing to accept sub-optimal solutions. // For a detailed algorithm see also "Intelligence artificielle" // (ia.txt, in French). For further local search strategies, // see "Softcomputing" (sc.txt, in German). Simulated Annealing for Graph Coloring: The application of simulated annealing to the graph coloring problem using the following formalization: * Be the set of colors {1,...K} * Be the graph (V,E) * The set of feasible solutions is the set of all mappings from nodes to colors, called "colorings" x:Node->{1,...K} * One defines for each color i=1..K the number of nodes getting this color under a coloring x: C[i](x) := | { v | v in V, x(v)=i } | * One defines for each color i=1..K the number of edges connecting two nodes of the same color under a coloring x: I[i](x) := | E /\ (C[i](x) x C[i](x)) | * The function to be minimized is the "penalty function" f(x) := 2 * SUM { C[i](x)*I[i](x) | i=1..K } - SUM { C[i](x)^2 | i=1..K } The first term of this function punishes illegal colorings and the second one favors large color clsses, i.e. few colors. * The neighborhood of a coloring is a coloring, which alters one color. Every local minimum of the function f is a feasible coloring. The idea used in this approach is that the set of feasible solutions is relaxed, such that it also contains illegal colorings. In turn, the penalty function f punishes illegal colorings so hard, that they will never be chosen. // This is the longest file I ever wrote. // May it help somebody out there.