Sunday, December 23, 2012

C/C++ code - performance optimization tricks

This post simply borrows (read cut-pastes) from material in
http://stackoverflow.com/questions/2074099/coding-practices-which-enable-the-compiler-optimizer-to-make-a-faster-program

http://andr0meda.com/plunge/?p=95

http://assemblyrequired.crashworks.org/2008/07/08/load-hit-stores-and-the-__restrict-keyword

They are my notes all in one place.I have not yet tried any of the code snippets myself. I take them at face value, to be tried as and when I'm called upon to do so.

 Memory Pools

Algorithms used to manage the heap are notoriously slow in some platforms (e.g., vxworks). Even worse, the time that it takes to return from a call to malloc is highly dependent on the current state of the heap. Therefore, any function that calls malloc is going to take a performance hit that cannot be easily accounted for. That performance hit may be minimal if the heap is still clean but after that device runs for a while the heap can become fragmented. The calls are going to take longer and you cannot easily calculate how performance will degrade over time. You cannot really produce a worse case estimate. The optimizer cannot provide you with any help in this case either. To make matters even worse, if the heap becomes too heavily fragmented, the calls will start failing altogether. The solution is to use memory pools (e.g., glib slices ) instead of the heap. The allocation calls are going to be much faster and deterministic if you do it right.

 For C++, may be writing or using custom allocators that employ pooling will work well.

Avoid Aliasing related slowdowns

Write to local variables and not output arguments! This can be a huge help for getting around aliasing slowdowns. For example, if your code looks like

Bad 
     void DoSomething(const Foo& foo1, const Foo* foo2, int numFoo, Foo& barOut)       
     { 
       for (int i=0; i<numFoo, i++)  
        { barOut.munge(foo1, foo2[i]); 
        }
     } 

the compiler doesn't know that foo1 != barOut, and thus has to reload foo1 each time through the loop. It also can't read foo2[i] until the write to barOut is finished. You could start messing around with restrict'ed pointers, but it's just as effective (and much clearer) to do this:

Good

    void DoSomethingFaster(const Foo& foo1, const Foo* foo2, int numFoo, Foo& barOut)
    {
       Foo barTemp = barOut;
       for (int i=0; i<numFoo, i++)
       {
         barTemp.munge(foo1, foo2[i]);
       }
       barOut = barTemp;
    } 
It sounds silly, but the compiler can be much smarter dealing with the local variable, since it can't possible overlap in memory with any of the arguments. This can help you avoid the dreaded load-hit-store (mentioned by Francis Boivin in this thread).

 

Load Hit Store issues

- read a from memory A into register X
- do something with a in X -> becomes b
- write b from register X to memory B
- read c from memory C into register Y 

There is nothing wrong with the code, nor with the reasoning behind it. But the CPU does have a problem when it has to execute it. The reason is that C and B can be equal in value, pointing effectively to the same address in memory. This means that c is in fact intended to be b (or a part of b – instructions sometimes chop binary values into lower and upper word fragments).
This dependency is recognized in the L1 cache, since all write instructions to memory temporarily block the CPU from reading from that cache until the value is written. So before the CPU can start reading from ‘C’, it has to wait until b is written to ‘B’, which could be at the very same memory address. [If the CPU would not wait, the value used as c is potentially not the one written by the previous instruction to memory at address B, and any further program flow / calculation might be off.] This waiting can take up quite some time, from 40 to 80 cycles, depending on your CPU. Despite advances in bus-logic and increased Front-Side-Bus performance, writing into memory is usually one of the slowest operations in the instruction set.
 
int CauseLHS(int *ptrA)
{
   int a, b;
   int *ptrB = ptrA; // B and A point to the same address
   *ptrA = 5; // write data to address ptrA
   b = *ptrB; // read that data back in again (won't be
                  // available for 80 cycles)
   a = b + 10; // stalls! the data isn't available yet.
} 

 

Use __restrict keyword for pointers

Use pointer parameters only if necessary, and mark them restrict if possible. This helps alias analysis a lot because the programmer guarantees there is no alias (the interprocedural alias analysis is usually very primitive). Very small struct objects should be passed by value, not by reference.
int slow(int *a, int *b)
{
   *a = 5;
   *b = 7;
   return *a + *b; // LHS stall: the compiler doesn't
                  // know whether a == b, so it has to
                  //  reload both before the add
}
 
int fast(int * __restrict a, int * __restrict b)
{
   *a = 5;
   *b = 7; // RESTRICT promises that a != b
   return *a + *b; // no stall; the compiler hangs onto
                         // 5 and 7 in the registers.
}
  
It bears repeating that __restrict is a promise you make to your compiler. If you break your promise, you can get incorrect results. If pointers pA and pB are __restrict, and pA == pB, that will cause mysterious bugs.
For example, in the fast() function above, if a does equal b, then the correct return value would be 14. However, the compiler won’t know that *b = 7 has changed the value of *a, and so you might actually get a return value of 12.

 

Pay attention to order in which memory is accessed

#define N 1000000;
int matrix[N][N] = { ... };

//awesomely fast
long sum = 0;
for(int i = 0; i < N; i++){
  for(int j = 0; j < N; j++){
    sum += matrix[i][j];
  }
}

//painfully slow
long sum = 0;
for(int i = 0; i < N; i++){
  for(int j = 0; j < N; j++){
    sum += matrix[j][i];
  }
}
  

 

Do IO in blocks

For example :- 
static const char  Menu_Text[] = "\n"
    "1) Print\n"
    "2) Insert new customer\n"
    "3) Destroy\n"
    "4) Launch Nasal Demons\n"
    "Enter selection:  ";
static const size_t Menu_Text_Length = sizeof(Menu_Text) - sizeof('\0');
//...
std::cout.write(Menu_Text, Menu_Text_Length);
 

 

Don't use printf family for constant data

Constant data can be output using a block write. Formatted write will waste time scanning the text for formatting characters or processing formatting commands. See above code example.

 

Declare small functions as inline or macros


 

 Loop Unrolling

Each loop has an overhead of incrementing and termination checking. To get an estimate of the performance factor, count the number of instructions in the overhead (minimum 3: increment, check, goto start of loop) and divide by the number of statements inside the loop. The lower the number the better.

unsigned int sum = 0;
for (size_t i; i < BYTES_TO_CHECKSUM; ++i)
{
    sum += *buffer++;
}
After unrolling:
unsigned int sum = 0;
size_t i = 0;
**const size_t STATEMENTS_PER_LOOP = 8;**
for (i = 0; i < BYTES_TO_CHECKSUM; **i = i / STATEMENTS_PER_LOOP**)
{
    sum += *buffer++; // 1
    sum += *buffer++; // 2
    sum += *buffer++; // 3
    sum += *buffer++; // 4
    sum += *buffer++; // 5
    sum += *buffer++; // 6
    sum += *buffer++; // 7
    sum += *buffer++; // 8
}
// Handle the remainder:
for (; i < BYTES_TO_CHECKSUM; ++i)
{
    sum += *buffer++;
}
In this advantage, a secondary benefit is gained: more statements are executed before the processor has to reload the instruction cache.

 

Use boolean arithmetic instead of if statements. 

Basically reduce if blocks.
For example:-

bool status = true;
status = status && /* first test */;
status = status && /* second test */;
 

 

Factor out variable allocation outside of loops


 

Declare static text (const literals) as static const

When variables are declared without the static, some compilers may allocate space on the stack and copy the data from ROM. These are two unnecessary operations. This can be fixed by using the static prefix.

 

Avoid too many floating point to int conversions and vice versa. 


They use different registers and converting from one type to another means calling the actual conversion instruction, writing the value to memory and reading it back in the proper register set.


 

Use const correctness as much as possible in your code. 

It allows the compiler to optimize much better.
In this document are loads of other optimization tips: CPP optimizations (a bit old document though)

Highlights:
  • use constructor initialization lists
  • use prefix operators (++i instead of i++)
  • use explicit constructors
  • inline functions
  • avoid temporary objects
  • be aware of the cost of virtual functions
  • return objects via reference parameters
  • consider per class allocation
  • consider stl container allocators
  • the 'empty member' optimization
  • etc

 

Put functions that are used together in close physical proximity in source.

 This increases the likelihood that they will be near each other in object files and near each other in your executable. This improved locality of instructions can help avoid instruction cache misses while running.

 

Make most functions static.

 

If global variables are used, mark them static and constant if possible.

If they are initialized once (read-only), it's better to use an initializer list like static const int VAL[] = {1,2,3,4}, otherwise the compiler might not discover that the variables are actually initialized constants and will fail to replace loads from the variable with the constants.

 

When writing tests with multiple conditions place the most likely one first.

 if(a || b || c) should be if(b || a || c) if b is more likely to be true than the others. Compilers usually don't know anything about the possible values of the conditions and which branches are taken more (they could be known by using profile information, but few programmers use it).

 

Using a switch is faster than doing a test like if(a || b || ... || z). 

Check first if your compiler does this automatically, some do and it's more readable to have the if though.

 

Data structure alignment

It  is the way data is arranged and accessed in computer memory. It consists of two separate but related issues: data alignment and data structure padding. When a modern computer reads from or writes to a memory address, it will do this in word sized chunks (e.g. 4 byte chunks on a 32-bit system). Data alignment means putting the data at a memory offset equal to some multiple of the word size, which increases the system's performance due to the way the CPU handles memory. To align the data, it may be necessary to insert some meaningless bytes between the end of the last data structure and the start of the next, which is data structure padding.

 

Profile, Profile, Profile

Then, and only then, you might want to start investigating the effects of telling the compiler how to use memory. Make 1 change at a time and measure its impact. Focus only on the hotspots.
Yes, compiler writers are not from a race of coding giants and compilers contain mistakes and what should, according to the manual and according to compiler theory, make things faster sometimes makes things slower. That's why you have to take one step at a time and measure before- and after-tweak performance.

And yes, ultimately, you might be faced with a combinatorial explosion of compiler flags so you need to have a script or two to run make with various compiler flags, queue the jobs on the large cluster and gather the run time statistics. If it's just you and Visual Studio on a PC you will run out of interest long before you have tried enough combinations of enough compiler flags.






 
 

Sunday, December 2, 2012

Testing the optimality of algorithms using Alloy - Part 3

Continuing from part 2...
http://feels-good-to-ramble.blogspot.in/2012/11/testing-optimality-of-algorithms-using_30.html
 , we are not done yet concluding that a BFS-based coloring algorithm is suboptimal. We now need to modify the algorithm to run BFS for every vertex as starting vertex in turn, and then check whether the least number of colors used in all those sub-runs is still suboptimal. This complicates things quite a lot, because we now need to store results for every BFS sub-run. It is not now enough to continue to store with Node, the fields that we stored with it in part 2. This is because, the field values are liable to change for each sub-run of the BFS. We solve this problem in a simple way by introducing another signature called Run, to denote a sub-run. We store fields, 'colormap', 'ordermap', 'levelmap' and 'snapshot', to denote maps between a Node and respectively its BFS 'color', its order of discovery during this sub-run of the BFS, its level in this sub-run of the BFS and finally, the  set of Nodes that include itself and the ones discovered before it in this sub-run. We need to introduce one more signature called Order which is the order of discovery of a Node for a sub-run. Only, the fields 'adjacency' and 'optcolor' remain with the Node, since they do not change from sub-run to sub-run. Finally, the predicate 'bfs_suboptimal' changes to assert the following:
all r:Run |  bfs[r] and #Node.optcolor < #r.colormap[Node]

which means that for every sub-run, the coloring found by Alloy (in 'optcolor') uses fewer colors than the coloring obtained by BFS in that sub-run.
  With these changes to our Alloy source,  it was found that no suboptimal coloring is obtained for vertices (Nodes)  up to 5 in number. The smallest graphs where BFS-based coloring is suboptimal has 6 or more vertices. All the test results in the figures below show only the 'optcolor' arrows and not the colors obtained through BFS. Also, the double-arrowed, red edges between any two Nodes (blue boxes) are the adjacency edges between the nodes.


Fig 1
Fig 2
Fig 3

 Some more results found after tweaking the source to constrain the number of 'optcolor's required to more than 3, then we get no 'suboptimal' coloring for 6 nodes, but we do get some for 7 nodes, one of which is shown in fig 4.

Fig 4

 This means that for 6 nodes the BFS based coloring does not use more than 4 colors in any case.




Fig 5

 It was also found similarly, that for 7 nodes the BFS based coloring does not exceed 5 colors, but for 8 nodes it does, as shown in Fig 5.
 The Alloy source listing is given below :


open util/graph[Node] as GR
open util/ordering[Order] as OR
open util/ordering[Color] as CO
open util/ordering[Level] as LV


sig Color {}

sig Run {
   start_node: one Node,
  private colormap: Node->one Color,
  private ordermap: Node one ->one Order,
  private levelmap: Node->one Level,
  private snapshot: Node  -> set Node
}

sig Order{}

sig Node{
    adjacency: set Node,
    optcolor: one Color
}

sig Level {}

pred  undirected_graph() {
    GR/undirected[adjacency]
    GR/noSelfLoops[adjacency]
       GR/stronglyConnected[adjacency]
}

pred init (r:Run) {
    let i=OR/first | let n= r.start_node | r.ordermap[n]=i and
    r.levelmap[n]=LV/first and r.snapshot[n] = n and r.colormap[n]=CO/first
}

fun done_nodes(r:Run, n:Node): set Node {
   r.snapshot[n]
}

fun next_level (r:Run, n:Node) : Level {
   LV/next[r.levelmap[n]]
}

fun prev_level (r:Run, n:Node) : Level {
   LV/prev[r.levelmap[n]]
}

fun level(r:Run, n:Node):Level {
   r.levelmap[n]
}

fun next_order(r:Run, n:Node): Order {
  OR/next[r.ordermap[n]]
}

fun prev_order(r:Run, n:Node): Order {
  OR/prev[r.ordermap[n]]
}

fun order(r:Run, n:Node) :Order {
  r.ordermap[n]
}

fun nodes_at_same_level(r:Run, n:Node) : set Node {
  (r.levelmap). (r::level[n])
}

fun node_at_order (r:Run, o:Order) : one Node {
  (r.ordermap). o
}

fun nodes_at_level (r:Run, l:Level) : set Node {
  (r.levelmap). l
}

fun next_siblings(r:Run, n:Node): set Node
{
  let parent_nodes= r::nodes_at_level [r::prev_level[n] ]   |
     some parent_nodes implies
      {
       let siblings = parent_nodes.adjacency -  r::done_nodes[n] |
         some siblings  implies
            siblings
         else
            none
     }
      else
            none
}


fun  children_at_next_level(r:Run, last_upper_node:Node) : set Node
{
    let n=last_upper_node |
           let peers= (r::nodes_at_same_level[n].adjacency - r::done_nodes[n]) |
           some peers implies
             peers
           else
             none
}

pred bfs(r:Run) {
    init[r]  
    all ord:Order - OR/last | step[r, ord ]
   correctness [r]
}

pred correctness [r:Run] {
   all n:Node | r.colormap[n] not in r.colormap[n.adjacency]
}

fact {
  some Node
  undirected_graph 
  no disj r1,r2:Run | r1.start_node=r2.start_node
  no disj n1,n2:Node | start_node.n1 = start_node.n2
  Run.start_node = Node
  coloring_prob
}

pred coloring_prob  { (no n,m:Node | (m in n.adjacency and n.optcolor= m.optcolor) ) }

pred set_color(r:Run, n:Node)
{
 let prevn = r::node_at_order[r::prev_order[n]] | let cols=Color - r.colormap[(n.adjacency & r::done_nodes[prevn])] |
   r.colormap[n] = CO/min[cols]
}

--This creates the next level from the current level. It is based on traversing the adjacencies
pred step (r:Run, i : Order) {
     let n= r::node_at_order[i] | 
     {
       let sibs=r::next_siblings[n] |
        some sibs implies
        {
         one m:sibs | r::order[m] = r:: next_order[n]  and
              r::level[m] = r::level[n] and
              r::done_nodes[m] = r::done_nodes[n] + m
       }
      else
      {
         let  children=r::children_at_next_level[n] |
         some children implies
         {
           one m: children |
               r::order[m] = r::next_order[ n ] and
              r::level[m] = r::next_level[n] and
              r::done_nodes[m] = r::done_nodes[n] + m
        }
      }
     let m= r::node_at_order[ r::next_order[n] ] |  r::set_color[m]
   }
}

pred bfs_suboptimal()
{
   all r:Run |  bfs[r] and #Node.optcolor < #r.colormap[Node]
}

run bfs_suboptimal for 6 but  4 Color

To conclude, in this set of 3 posts we saw how Alloy can be employed to test our understanding of problems, especially when we are learning about or exploring a subject. Specifically, we demonstrated how a simple (class-P) BFS-based algorithm is unable give optimal results for the vertex coloring problem, for many graph instances, shattering any intuition that this problem may have this algorithm as a solution. Of course, it does not demonstrate that there is no exact solution for the vertex-coloring problem in P, simply that our BFS-based algorithm is not the one.
 Of course, the world knows the vertex-coloring problem is NP-complete, and that there is perhaps (note the 'perhaps') little sense in pursuing an exact solution in P for this problem. These posts simply illustrate a larger point about how we can use approaches and tools to test the efficacy of any approximate algorithms that we may rig up for such problems. For example in future it may be interesting to find out the bounds on the suboptimality of the BFS-based algorithm, that is, how far from the optimal solution can it get, and whether there are bounds on this distance from optimality. This requires theoretical proof, but tools like Alloy can suggest directions and check hunches.

Friday, November 30, 2012

Testing the optimality of algorithms using Alloy - Part 1

The Sea from Candolim beach, Goa
  I'm just back after a wonderful holiday in sunny Goa. Before embarking on that trip, I was nosing around some topics related to Theory of Computation, Graphs and Graph coloring. Sound like ponderous names ? Unfortunately, they are at the heart of some very difficult yet practical problems that experts try their best to solve satisfactorily and with a reasonable amount of computer resources (time and memory). They are what are called intractable problems, as against problems (called class-P  or polynomial time complexity problems) like solving a set of linear equations using Gaussian elimination or finding the shortest path between two cities, given a network of roads between cities and the distance between any two cities. These intractable problems span the entire gamut of disciplines that have contributed to all aspects of modern living - job scheduling problems, traveling salesman problems, optimal register allocation during compilation (so that the compiler can generate code that runs fast), 0-1 integer programming, Set covering etc. See the paper by Richard Karp that mentions  several more of them - http://www.cs.berkeley.edu/~luca/cs172/karp.pdf



Fig 1

   I was also playing with a model checking tool called Alloy (http://alloy.mit.edu/alloy/). What exactly is Alloy ? Let me quote from the Alloy site...
"Alloy is a language for describing structures and a tool for exploring them. It has been used in a wide range of applications from finding holes in security mechanisms to designing telephone switching networks.
An Alloy model is a collection of constraints that describes (implicitly) a set of structures, for example: all the possible security configurations of a web application, or all the possible topologies of a switching network. Alloy’s tool, the Alloy Analyzer, is a solver that takes the constraints of a model and finds structures that satisfy them. It can be used both to explore the model by generating sample structures, and to check properties of the model by generating counterexamples. Structures are displayed graphically, and their appearance can be customized for the domain at hand.
...

The Alloy Analyzer works by reduction to SAT. "

 (SAT means the "Boolean Satisfiability Problem", which, briefly, means determining whether or not there exist true/false assignments to boolean variables that satisfy given constraints or boolean formulas made from 'and' and 'or' boolean, binary operations between those variables. See http://en.wikipedia.org/wiki/Boolean_satisfiability_problem. For the purpose of this post, the implication of the last line of the quoted paragraph above, is that the model of any intractable problem, that we have built in the Alloy language, can be analyzed , with some limitations that will be clarified in later posts, by the Alloy analyzer after first converting (reducing) the model to boolean formulas and applying available analyses for SAT problems. To reinforce the words of caution, I am not saying that Alloy can solve a problem for any input.)

I had used Alloy in some work earlier, for modeling CAD to PLM data conversion situations, the purpose being to discover any inconsistencies and discrepancies that may arise because of the way we modeled. These gotchas  sometimes ( or is that 'often' ?) escape us when we are trying to use only our grey cells to analyze and check complex models to see if they behave the way we expected or  to see whether there are instances that behave the way we don't want them to.


Graph coloring

  This time my requirement was slightly indulgent. During the course of some reading on Graph algorithms, I came across the line that 'The Graph coloring problem is NP-complete' and I was puzzled why. An oft-cited application of graph coloring is the problem of optimal register allocation for local variables and function arguments for a function call, during the code generation phase of compilation. The compiler ( its code generation/optimization phase) will first create a graph whose vertices correspond to the variables, and two vertices have an edge between them if their lifetimes overlap, which implies that the two variables can't be allocated the same register during code generation. It is important to know the minimum number of registers needed, because their total number is limited, and we would like as many of the variables to be accommodated in them, since access to registers is fastest for the CPU. First, what is the graph coloring problem, and specifically vertex coloring ?
 This is what Wikipedia has to say.
In graph theory, graph coloring is a special case of graph labeling; it is an assignment of labels traditionally called "colors" to elements of a graph subject to certain constraints. In its simplest form, it is a way of coloring the vertices of a graph such that no two adjacent vertices share the same color; this is called a vertex coloring.


Fig 2



 

 Graph coloring and the classes NP, NP-hard and NP-complete

  One flavor of the problem is to find the least positive integer k, representing the number of colors, that can vertex-color a graph (such that no two adjacent vertices share the same color). This is an optimization problem  and this problem is said to be NP-hard. We define this term in the process of defining another class of problems, called the class of  NP-complete problems. In lay terms, the last term, NP-complete refers to a class of problems and means two things :-

 1) that any problem of this class belongs to the class of NP  (non-deterministic, polynomial complexity) problems
This in turn means that the problem has a known solution that will, in the worst case, take a number of steps (to completion) that is an exponential function of the problem size ( size being in this case the number of vertices of the graph, for example) on a deterministic RAM (Random access machine - pretty much means the kind of machines we use) and involves evaluating an exponentially large number of branches of execution. But the same problem can finish in polynomial time on a fictitious non-deterministic machine which can multi-task and work several branches of execution concurrently to get the solution. 'Exponential' means 'impractical to solve' on any Computer (however powerful) based on the RAM concept. In short, exponential time means 'BAD'.

 2) that any problem of this class belongs to the class NP-hard
This means that any problem in NP can be reduced 'easily' (i.e. in polynomial time) to this problem.Without  further delving into what this exactly means, we shall be content with what this implies - that if someone is clever enough to find a solution for even one  problem in this class (NP-hard), in polynomial time ( which means the algorithm becomes closer to being easily or practically implemented on computers; in short, it's a  'GOOD' algorithm), then all the other (yes, every other) NP problem can also be solved in 'GOOD' time too.Since the consequent of the earlier 'if' sentence is believed to be false by most experts, the antecedent ( that we'll find a polynomial time algorithm for an NP-complete problem some day) is also very likely to be false.

 In other words, an NP-complete problem is an NP problem that is also an NP-hard problem.
   In fact, even the seemingly simpler decision (yes/no) problem of checking whether k colors (k>2) are enough to vertex-color a graph, is claimed to be NP-complete.

Proving that Graph coloring is in NP 
 Well, at least proving that the decision problem is in NP is easy.  To do this, we need to demonstrate that there is an algorithm that solves it taking exponential time on a deterministic RAM, but polynomial time on a non-deterministic RAM. Let's devise such an algorithm as follows (we are talking about undirected graphs here, which means that if vertex a is connected to vertex b, it implicitly means that vertex b is also connected to vertex a).  Picture that we have k bins into which we want to distribute the n vertices of the graph. The distribution is such that all edges of the graph are only between 2 vertices in different bins. Going from an arbitrarily chosen first vertex to the nth, we can choose a bin in k ways for the first vertex, in k ways for the second and so on. So, this way we actually evaluate k^n (k to the n) vertex to bin assignments in the worst case. Evaluation of a candidate solution involves checking if any vertex's color assignment violates the edge constraint mentioned earlier. If it does, then the solution is rejected. If there is no violation for any vertex, then this is an acceptable solution and we can stop. In the worst case we may need to evaluate all possible solutions (k^n) before we can pronounce a 'no'. This no-brainer approach is called a 'brute-force' algorithm and its time complexity is proportional to k^n. This algorithm is said to be of order O(k^n), i.e. exponential complexity. Moreover, if we had a hypothetical machine that could work by branching out at each decision point of choosing a bin for a vertex, and all the branches could concurrently execute,  then such a non-deterministic machine could find an answer (yes/no) in polynomial time, this time being the time taken by a single branch, which is O(n). By that token, this problem belongs to the complexity class NP.
  Note that if k were 1, the coloring succeeds only if the graph is an n-vertex graph with n components (i.e. the graph has no edge at all). So the algorithm to see if the graph is 1-colorable is trivial - you just check if there is no edge at all. If  k=2 it means we need to check if the graph has a 2-coloring. The 2-coloring succeeds only if the graph is bipartite i.e. its vertices can be divided (partitioned) into two disjoint sets (bins) such that any edge of the graph is only between a vertex of one set and that of the other set ( and not between two vertices of the same set).
 Now, the 2-coloring problem ( or equivalently, the problem of determining whether a graph is bipartite) can easily be solved in polynomial time by making a breadth-first search on the graph (BFS), which divides the vertices into different levels that represent the edge distance from the starting vertex chosen for the search, and checking whether there is any edge between vertices at the same level. If there is an edge between two vertices at the same level during a BFS on the graph, then we can say that the graph is not bipartite and, equivalently, cannot be 2-colored. The 1-coloring and 2-coloring problems are said to belong to time complexity class P (for polynomial complexity). This is because the method used, BFS, has a time complexity of order O(m+n), m being the number of edges, and m+n is, in the worst case, a polynomial function in n ( the maximum possible value of m is order O(n^2), so m+n is at most order O(n^2)).
 For BFS level assignment, you select a starting vertex arbitrarily, assign it level 0, then assign all vertices adjacent to it , a level 1, and all vertices adjacent to the level 1 vertices, a level 2, and so on, till all unassigned vertices are assigned a level. It is a feature of the BFS algorithm, that the level of any vertex is either exactly one less, or one more, or the same value as any vertex adjacent to it. If and only if a graph is bipartite, can we assign some color A to the even level vertices and color B to the odd level vertices and have a valid 2-coloring.
 Why can we or can't we extend this algorithm to see if a graph is 3-colorable or higher number-colorable ? Suppose that during BFS coloring, if we find two connected vertices at the same level, we simply choose another color. For example, in the above Fig 2 showing BFS-based coloring, we see this situation occurring for vertex 3 and then vertex 4, both at level 1, hence requiring vertex 3 to be colored Blue instead of Green, and then vertex 4 being colored Black. 

Is the BFS-based approach optimal ?
So we get a nagging doubt here whether this BFS-based coloring process might result in a situation where we might wastefully need say, an extra color in later steps, because of the order in which we chose the starting and other nodes for assigning colors. Putting it differently, if we needed to use say, 4 colors using a BFS-based coloring approach, we wonder whether some other order would in fact have given us a better coloring e.g. a 2- or 3-coloring. On the other hand, may be some other policy for assigning a color to a vertex could have given a better solution i.e. a 2- or 3-coloring for example. In short, if, following the above BFS-based coloring algorithm, we did need a 4th color, then how can we be sure that there is no 2- or 3-coloring for this problem following some other algorithm, and that we indeed require at least 4 colors for vertex-coloring the graph irrespective of how we went about it ?
 We need to be able to prove that our algorithm ( the BFS-based one described above was just an example, probably a very simple-minded one) always gives the best solution for any graph; or we need to produce a graph as a counter-example to show that our algorithm is not good enough because there is a coloring for the graph that uses fewer colors. Although not impossible, it is non-trivial to prove either without some help. Here is where Alloy comes in. Let's look at that in my next post
 


Testing the optimality of algorithms using Alloy - Part 2

 Continuing where we left off last time on
http://feels-good-to-ramble.blogspot.com/2012/11/testing-optimality-of-algorithms-using.html
, I give the description of a first-shot Alloy model that tries to check whether there is a graph for which the Alloy analyzer can find a better (cheaper) vertex coloring than the one based on BFS (Breadth First Search).
It is not the intent of this post to be a tutorial for Alloy. There are more competent places for that viz.
  http://alloy.mit.edu/alloy/documentation.html
But some words to build the context are in order.  Pardon me if this is turning out to be some kind of 'notes to myself' article.
 How many times has it happened that we are confident that an algorithm we devised is correct, only to find some days later, and sometimes after the solution reaches the Testing department or the Customer, that there was one or more important scenarios that we missed out leading to costly recoding/fixing/rework not to mention embarrassment ?
 When the nature of the problem is such that there are a combinatorially large number of possible cases to test, Test engineers may also miss some or all of these cases, just because there are so many of them even for a relatively small problem size, and the tests simply did not cover the problem cases ( which might be like needles in a haystack). Now, consider the situation where these one off trouble cases happen to introduce say, a possible security threat, or cause the application to hang completely in a deadlock, or give a completely way-off calculation for the insurance premium to be charged. Then, it becomes important for Testing to cover all cases, so as to instil confidence in the provided solution. The Alloy tool helps in this humanly tedious (if not impossible) task, by an ability to deliver at the analysis and design (i.e. pre-coding) stage itself. Developers can model their solution at a high level of abstraction in the Alloy language, and run checks on what they have modeled, long before actual coding in the implementation language. Alloy lets you test the model in two ways.
  •   produces instances of the model or tells you the model is inconsistent and no instance exists. This is done via the 'run' command
  •  produce counter-examples for an assertion coded in Alloy. This is done via the 'check' command. If the assertion is a valid one then Alloy cannot produce any counter-examples and reports likewise.
Both the above tests are done for a limited scope, i.e. all the combinatorially many solutions are evaluated for only a fixed number of objects of a 'signature' i.e. set (the phrase in bold corresponds to the phrase, 'small problem size' above. In Alloy 'fixed' is not necessarily very small, because SAT analyzers today can address fairly large problems, even in the case where the problem is NP-complete). The 'run' command might be specified like this in the Alloy file.

run bfs_suboptimal for 10


where bfs_suboptimal is the name of a predicate that we defined earlier in the model. In this case, 10 is the fixed number of objects of any signature that the Alloy analyzer will process. The basis for trusting the analysis for a fixed number, or 'scope' of the analysis, is that most often the results carry over to higher numbers. But one needs to 'experiment' with several values before achieving some confidence level.

The 'check' command might be specified like this.

check coloringIsCorrect for 6

where coloringIsCorrect is the name of an assert block that was defined in the model.

First remember that the Alloy modeling language is a declarative language after the style of Prolog or Lisp. This means the statements of an Alloy program describe the constraints etc. that hold in the resultant model, and the Alloy analyzer derives the model that satisfy these constraints. This is a paradigm shift for most of us programmers who have been raised on a staple of imperative languages like Java, C, C++ etc. where the program statements are actually instructions about what steps to take, the resulting constraints and values etc. being only implied. To mention one big difference, Alloy has no 'assignment' statement. Something like "x=x+1" is simply a 'FALSE' statement in Alloy, just as in Mathematics, and does not assign anything to x. Only some familiar imperative languages have accorded legitimacy to such syntax. To give an example of an Alloy model ( we use the term 'model' interchangeably for an Alloy program source as well as for the structures generated from the program by the Alloy analyzer). Let's take a directed acyclic graph (DAG) data structure. This would be specified in Alloy as follows :

--Directed Acyclic Graph
sig Node{
   children: set Node
}

fact {
--The first fact asserts that there is at least one node in our model
  some Node
--The next fact asserts that no node is its own ancestor, i.e. that the structure
-- has no cycles
  no node:Node | node in node.^children
--The next fact asserts that any two nodes have a path from at least one of them
--to the other
  all disj n1,n2:Node | n1 in n2.^children or n2 in n1.^children
}

--Predicate to keep generating and showing instances of the data structure 
--described by the above model
pred show(){}

--Run the predicate for at most 5 nodes
run show for 5

  Some of the DAGs obtained by running the analyzer on the above model are shown in the images of figures 1a to 1d.
Fig 1a
Fig 1b
Fig 1c
Fig 1d
For our problem of deducing whether the BFS-based vertex-coloring algorithm we outlined in the last post is suboptimal, our approach will be to define a predicate called bfs_suboptimal which asserts that the number of colors as derived by the Alloy analyzer, given a free hand so to say, is less than the number of colors that are required if the coloring follows the BFS-based approach. The BFS-coloring is described in the predicate bfs. This is a bit tricky to code, because the BFS algorithm has more of a "do this-do that" nature, and mapping that to the declarative paradigm of Alloy can be non-trivial. We look at this algorithm as follows - note we use the terms 'vertex' and 'node' interchangeably. First we choose a starting vertex (Node), which will form the root of the BFS tree to be generated by the algorithm, and assign it the first color in the Color set. We then 'discover' nodes adjacent to this node and color them, and then the (as yet undiscovered) nodes adjacent to those in turn and so on. At any stage, we use the fact that all nodes adjacent to a given node, except for the ones that have already been discovered in a previous stage by the BFS algorithm,  are at the same distance from the starting vertex, i.e. at the same 'level'. Any node will be then given a color that is the 'lowest' color  in the set of colors that were earlier assigned to all discovered nodes that are adjacent to this node (for example if all colors are numbered 0,1,2,...,  and the set of colors contains the colors 2,4,8 then color 2 is the lowest in this set). Before we proceed further, I must clarify something here. The use of words 'earlier', 'discover' etc. gives a connotation of time, which can be misleading. As mentioned earlier there is no assignment, so it is not as if a single set of discovered nodes get updated from stage to stage as new nodes are discovered. On the contrary, we need to store (this is one way to do it) what may be called a snapshot ( set) of discovered nodes at every stage, unique to that stage. I'll  provide the details as we discuss the program itself.

We first include several modules that we will use, and that come with the Alloy installation. These are parametrized modules and they define predicates related to Graphs, and predicates related to there being a total ordering defined for Nodes, Colors and Levels.
open util/graph[Node] as GR
open util/ordering[Node] as trnode
open util/ordering[Color] as CO
open util/ordering[Level] as LV

Then we define several 'signatures' which are like Java Classes, and also like Sets of Mathematics. The signature 'Node' represents the set of vertices of a graph. 'Level' represents the set of levels. The level of a vertex is its least distance from the starting vertex when a BFS tree is created from the graph. The signature Node has some interesting fields. The field 'adjacency' is the set of Nodes adjacent to a Node. The field 'done_nodes' is the snapshot set of Nodes discovered at the end of a stage. Since at every BFS stage exactly one new Node is discovered, we store this snapshot with a Node. The field 'color' is a Node's color as per the BFS algorithm, while the field 'optcolor' is the color as per only the bare minimal constraints defined i.e. as per the non-BFS solution.

sig Color {}

sig Node{
    adjacency: set Node,
    level: Level,
    done_nodes: set Node, 
    color: one Color,
    optcolor: one Color
}

sig Level {}

 This is followed by several predicates, most of which are internal ones, but some of them are also later used as  predicates for the 'run' command.

This predicate asserts  (by the way, see http://www.thefreedictionary.com/predicate for meanings of the word 'predicate') that elements of the set Node form an undirected graph with no self loops and that the graph is strongly connected i.e. that all of the undirected graph is just one connected component.

pred  undirected_graph() {
    GR/undirected[adjacency]
    GR/noSelfLoops[adjacency]
       GR/stronglyConnected[adjacency]
}


This next one is an Alloy function (a function in Alloy simply wraps an Alloy expression). It returns  a sibling, if any,  at the same level as the input Node n - this requires you to look at the graph with a BFS tree superimposed on it, the starting vertex for the BFS being the root of the BFS tree. Since there may be one or no such sibling in the BFS tree, the return value of the function is qualified with the keyword 'lone', meaning 0 or 1.

fun next_sibling(n:Node): lone Node
{
  let parent_nodes= level.(LV/prev[n.level]) |
     some parent_nodes implies {
       let siblings = parent_nodes.adjacency -  n.done_nodes |
        some siblings  implies
          trnode/min[siblings]
        else
          none
     }
    else
      none
}

This Alloy function will (as is obvious from its name) find the child at the next level, from input Node, if any.

fun first_child_at_next_level(last_upper_node:Node) : lone Node
{
   trnode/min[last_upper_node.level.~level.adjacency - last_upper_node.done_nodes]
}

The following predicate is the workhorse for the BFS algorithm. It first calls the init predicate explained later. Then it calls the step predicate, passing two consecutive levels (remember Levels are totally ordered). This is very typical of how we specify in Alloy how two stages of an algorithmic model are related. Finally it imposes the correctness constraint for vertex coloring, so as to further narrow down the graphs enumerated by Alloy to only the valid ones.
pred bfs() {
    init 
    all n:Node - trnode/last  |    { let n'=trnode/next[n]  | step [n,n'] } 
    correctness
}

The following predicate initializes the BFS, by asserting that the first Node ( there is a total order on Node too) gets the first color ( as calculated by the BFS algorithm) from Color (there is a total order on Color too), and optcolor, which represents the color as calculated by Alloy (the prefix 'opt' stands for optimal, which is something of a misnomer), is also initialized to the same color.

pred init () {
    let n=trnode/first |
    n.level=LV/first and n.done_nodes=n and n.color=CO/first and n.optcolor=n.color
}

This predicate asserts the basic vertex-coloring requirement that any coloring produced by BFS must satisfy, i.e. that a Node's 'color' must be different from that of its adjacent Nodes.


pred correctness {
   all n:Node | n.color not in n.adjacency.color
}

This assert is just a wrapper over the predicate 'correctness', so that we can run some 'check' commands (as mentioned before) to ensure the generated models have correct coloring.

assert coloringIsCorrect {
    correctness
}

The following predicate asserts the basic vertex-coloring requirement that any coloring produced using the 'optcolor' field must satisfy, i.e. that a Node's 'optcolor' must be different from that of its adjacent Nodes. Note how this depends only on specifying the bare minimum constraint. This predicate is used to generate the 'optcolor' fields. Being as minimally constrained as it is, this predicate will not actually generate optimal colorings if used in isolation, but the single predicate 'bfs_suboptimal', described later will let us do what we set about to i.e. to prove that it generates a coloring that is better than the BFS-based, thus showing that BFS-coloring is suboptimal.


pred coloring_prob  { (no n,m:Node | (m in n.adjacency and n.optcolor= m.optcolor) ) }

The following predicate sets the BFS 'color' of a Node to be the minimum from the set of Colors that is the complement of the BFS 'color's of the already discovered Nodes adjacent to it - that's a mouthful, but you get the idea. The effect is that the input Node is given the least color chosen from the entire set of Colors, minus the ones given to adjacent Nodes that have been discovered thus far.

pred set_color(n:Node)
{
 let cols=Color - (n.adjacency & n.done_nodes).color |
   n.color = CO/min[cols]
}

The following predicate is probably the most important one for the BFS coloring. It relates two consecutive BFS coloring stages. Since a stage corresponds to a Node, it takes two Nodes  as arguments. These Nodes are the ones that were discovered consecutively in the two consecutive stages. The second argument Node is asserted to be either a sibling (which here means either an immediate sibling or a cousin) of the first argument Node or the first child at the next level. Finally, the second Node is colored using the predicate 'set_color' described earlier.

--This creates the next level from the current level. It is based on traversing the adjacencies
pred step (n,n': Node) {
     some n::next_sibling implies
     {
       n'=n::next_sibling and n'.level= n.level and n'.done_nodes = n.done_nodes + n'
     }
    else {
     let m=n::first_child_at_next_level  |
      some m implies {
       n'=m and n'.level = LV/next[n.level] and  n'.done_nodes = n.done_nodes + n'
      }
      else {
         n'.level = n.level
         n'.done_nodes = n.done_nodes
      }
    }
    set_color[n']
}

This fact brings it all together, ensuring that the Graph structures generate by Alloy have Node 'color's as per BFS, and Node 'optcolor's as per the bare minimum constraints as specified in predicate 'coloring_prob'.
fact {
  undirected_graph 
  bfs
  coloring_prob
}
This predicate further constrains our model to those in which the number of colors used by the BFS is higher than the ones obtained from 'coloring_prob'. 
If there happens to be some Graph structure that Alloy finds for this predicate, then it means our BFS-based  coloring algorithm is suboptimal for those graphs, and therefore suboptimal in general.

pred bfs_suboptimal()
{
   #Node.optcolor < #Node.color
}
 Finally, we run the predicate 'bfs_suboptimal' for at most 5 objects of any signature, except for the signature Color, which we limit to 4. I actually started with values 3,then 4 etc. for the generic scope value, but didn't find any suboptimal graph coloring. The first suboptimal graphs started appearing for values 5 onwards. Similar experimentation was conducted for the scope value of Color.

run bfs_suboptimal for 5 but 4 Color

The following images are for some of the instances produced by the Alloy analyzer for the above 'run'. The red, double-arrowed edges between the blue Nodes denote the adjacency. A saffron colored arrow points from a Node to the 'optcolor' set for a Node, while any brown arrow points from a Node to the BFS 'color' for that Node. In all the instances, Node0 is the starting vertex for the BFS.

Fig 2a

Fig 2b

Fig 2c
Here's a fairly complex graph of 10 Nodes for a different run. I have eliminated the Color objects to avoid clutter.

Fig 2d
 So we seem to have found a fairly dependable way to check whether our algorithms behave as expected (e.g. whether our BFS-based coloring is suboptimal, in this case). But are you sure ? The alert reader may argue that our model is flawed. Our BFS algorithm was run only for a single Node as starting vertex. May be if we had coded our Alloy program to run the BFS for every node of a graph as starting vertex in turn, and then called 'bfs_suboptimal', may be we would never have found any suboptimal graph. Note that this way our algorithm would still be in complexity class P, since the time would simply get multiplied by n, the number of vertices. 
 We'll see about this in the next post.

Is this stock advise worth taking seriously?

 Introduction Business related TV channels and newspapers are replete with stock advisory from investment firms and certified individua...