Datacenter Cooling (or Counting Hamiltonian Paths in a Grid)

// Datacenter Cooling (or Counting Hamiltonian Paths in a Grid)
//
// Written by Grant Jenks
// http://www.grantjenks.com/
//
// DISCLAIMER
// THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE
// LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
// OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND,
// EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
// ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.
// SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY
// SERVICING, REPAIR OR CORRECTION.
//
// Copyright 2012
// This work is licensed under the
// Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License
// To view a copy of this license, visit
//    http://creativecommons.org/licenses/by-nc-nd/3.0/
// or send a letter to
//    Creative Commons
//    171 Second Street, Suite 300
//    San Francisco, California, 94105, USA
//
// I originally saw this problem at http://www.quora.com/challenges, but it
// looks like an ACM competition / AI textbook problem.
//
// Build with:
//    cl.exe /EHsc /O2 [/GL] [/DMEASURETIME] datacenter.cpp
//
//                             Datacenter Cooling
// We have some rooms in our datacenter, and we need to connect them all with
// a single cooling duct.
// Here are the rules:
//     * The datacenter is represented by a 2D grid.
//     * Rooms we own are represented by a 0.
//     * Rooms we do not own are represented by a 1.
//     * The duct has to start at the air intake valve, which is represented
//       by a 2.
//     * The duct has to end at the air conditioner, which is represented
//       by a 3.
//     * The duct cannot go in multiple directions out of the intake or the
//       AC - they must be the two endpoints of the duct.
//     * The duct must pass through each room exactly once.
//     * The duct cannot pass through rooms we do not own.
//     * The duct can connect between rooms horizontally or vertically but
//       not diagonally.
// Here is an example datacenter:
//     2 0 0 0
//     0 0 0 0
//     0 0 3 1
// There are two possible ways to run the duct here:
//     2--0--0--0
//              |
//     0--0--0--0
//     |
//     0--0--3  1
// or
//     2  0--0--0
//     |  |     |
//     0  0  0--0
//     |  |  |
//     0--0  3  1
// Write a program to compute the number of possible ways to run the duct. For
// the above example, the correct answer is 2.
//
// Input format:
// Your program should read from stdin. The input will be a series of
// integers separated by whitespace.
// The first two integers will be W and H, with width and height of the
// datacenter.
// These will be followed by W*H more integers, specifying the 2D grid
// representing the datacenter.
//
// Output format:
// Your program should write a single integer out to stdout: the number of
// possible ways to run the duct.
//
//                                 Test Inputs
// Small (result: 2)
//    4 3
//    2 0 0 0
//    0 0 0 0
//    0 0 3 1
// Medium (result: 20)
//    5 4
//    2 0 0 0 0
//    0 0 0 0 0
//    0 0 0 0 0
//    0 0 0 0 3
// Benchmarking (result: 301716)
//    7 8
//    2 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    0 0 0 0 0 0 0
//    3 0 0 0 0 1 1
//
//                                  Changelog
// 2012-03-15 Grant Jenks
// * Did a little research online and learned a new trick. If, while computing
//   reachability (the inner-dfs loop), a dead-end is discovered and the
//   dead-end is not the end position, then there are no possible paths. This
//   is a cheaper and narrower test for the observation in 2010-07-24.
//   Processing time: ~2.2 seconds
//
// 2012-03-14 Grant Jenks
// * Realized that I could memoize the tape. What if you occupy the same rooms
//   in two different configurations? Wouldn't the result be the same after
//   that point? For example:
//   If the starting board is:
//   S 0 0 0 0
//   0 0 0 0 0
//   0 0 0 0 0
//   0 0 0 0 0
//   0 0 0 0 E
//   And we reach this point (numbers indicate traversal pattern):
//   S 1 8 0 0
//   3 2 7 0 0
//   4 5 6 0 0
//   0 0 0 0 0
//   0 0 0 0 E
//   Then we will compute how many paths are reachable.
//   And if later we reach this point:
//   S 7 8 0 0
//   1 6 5 0 0
//   2 3 4 0 0
//   0 0 0 0 0
//   0 0 0 0 E
//   Then we know already (from memoization of the results) that this will
//   yield the same number of paths.
//   I prototyped these changes in a Python version of this program. While
//   the memoized results were frequently hit, cache maintenance and size
//   were impractical.
//
// 2010-07-24 Grant Jenks
// * Came up with an idea to approximate this problem as a max-flow problem.
//   This led me to consider what would happen if there existed an open
//   room, r, not adjacent to the end room, e, such that closing the room
//   would make the DFS test fail. Before this is done, I think it would
//   be good to add more comments and refactor the code some.
//   I got distracted from completing the above and didn't come back to the
//   code until a different idea come to me two years later.
//
// 2010-04-03 Grant Jenks
// * Converted conditionals in loop to straight line code.
//   Processing time: ~40 seconds
//
// 2010-03-27 Grant Jenks
// * Added BOUNDARY to datacenter array to eliminate bounds checking.
//   Processing time: ~60 seconds
//
// 2010-03-20 Grant Jenks
// * Added post-inner-DFS loop check that all open rooms were reachable.
//   Processing time: ~2 minutes
//
// 2010-03-13 Grant Jenks
// * Added the inner DFS loop but only checked that end was reachable.
//   Processing time: ~2 hours
//
// 2010-03-06 Grant Jenks
// * Created the basic outer loop with tape.
//   Processing time: Several hours.
 
#define OPEN 1
#define CLOSED 0
#define VISITED 0
#define BOUNDARY 0
#define DFS_VISIT 8
#define VISIT_MASK 0x80000000
 
// For data input/output
 
#include <iostream>
#include <iterator>
 
// For measuring time accurately
 
#ifdef MEASURETIME
#include <windows.h>
#endif MEASURETIME
 
// To coerce the compiler to use certain forms
 
#include <intrin.h>
 
// Unsigned int sizes
 
typedef unsigned __int8 __uint8;
typedef unsigned __int16 __uint16;
typedef unsigned __int32 __uint32;
 
using namespace std;
 
int main()
{
#ifdef MEASURETIME
  LARGE_INTEGER start_time;
  QueryPerformanceCounter(&start_time);
#endif
 
  __uint32 o_width, o_height, a_width, a_height;
  __uint32 datacenter_size, start, end, not_owned;
 
  // Use __uint16 though the backing array will be __uint8 so that input is not
  // interpreted as characters.
 
  typedef istream_iterator<__uint16> iis;
 
  // Read in the initial width and height and initialize the datacenter array.
 
  iis input = iis(std::cin);
 
  o_width = *input++;
  o_height = *input++;
  a_width = o_width + 2;
  a_height = o_height + 2;
  datacenter_size = a_width * a_height;
  __uint8 * datacenter = new __uint8[datacenter_size];
 
  // The perimeter of the datacenter is denoted by the BOUNDARY. This allows us
  // to skip bounds checking when looking at adjacent rooms to visit.
 
  for (__uint32 pos = 0; pos < a_width; pos++)
    {
      datacenter[pos] = BOUNDARY;
      datacenter[pos + (a_width * (a_height - 1))] = BOUNDARY;
    }
  for (__uint32 pos = 0; pos < a_height; pos++)
    {
      datacenter[pos * a_width] = BOUNDARY;
      datacenter[(pos * a_width) + (a_width - 1)] = BOUNDARY;
    }
 
  // Read in the rest of the datacenter.
 
  start = 0;
  end = 0;
  not_owned = 0;
 
  for(__uint32 o_pos = 0; input != iis(); input++, o_pos++)
    {
      __uint8 value = *input;
      __uint32 o_pos_w = o_pos % o_width;
      __uint32 o_pos_h = o_pos / o_width;
      __uint32 pos = (o_pos_w + 1) + ((o_pos_h + 1) * a_width);
 
      switch (value)
	{
	case 0:
	  datacenter[pos] = OPEN;
	  break;
	case 1:
	  not_owned++;
	  datacenter[pos] = CLOSED;
	  break;
	case 2:
	  if (start != 0)
	    {
	      cerr<<"Start defined twice!"<<endl;
              return 1;
	    }
	  start = pos;
	  datacenter[pos] = OPEN;
	  break;
	case 3:
	  if (end != 0)
	    {
	      cerr<<"End defined twice!"<<endl;
	      return 2;
	    }
	  end = pos;
	  datacenter[pos] = OPEN;
	  break;
	default:
	  cerr<<"Invalid input: "<<value<<endl;
	  return 3;
	}
    }
 
  // Define the tape.
 
  // Note that we're playing a game with the visited bit by putting
  // it in the high-bit of the tape position. The positions stored
  // on the tape are most likely less than 2^31-1 so this won't be
  // an issue.
 
  __uint32 * tape = new __uint32[4 * datacenter_size];
  __uint32 * tape_pos = tape;
  __uint32 * dfs = new __uint32[datacenter_size + 1];
 
  // Because we've recorded the 'start' and 'end' positions, we can
  // open the rooms.
 
  datacenter[start] = OPEN;
  datacenter[end] = OPEN;
 
  *tape_pos = 0x7FFFFFFF;
 
  // Rather than recursing, the tape is used to solve the problem
  // iteratively. This gives finer-grain control over the assembly.
 
  tape_pos++;
  *tape_pos = start;
 
  __uint32 steps = 0;
  __uint32 routes = 0;
  __uint32 max_steps = o_width * o_height - not_owned;
 
  while (tape_pos != tape)
    {
      __uint32 pos = *tape_pos;
 
      if (pos == end)
	{
          steps++;
	  *tape_pos |= VISIT_MASK;
	  if (steps == max_steps)
	    {
	      routes++;
	    }
	}
      else
	{
	  __uint32 * stack_pos = tape_pos;
	  __uint32 * dfs_pos = dfs;
 
	  int visited_count = 0;
 
          // Do a DFS through the open rooms to make sure all are
          // reachable and none (but the end room) are a dead-end.
 
	  while (stack_pos >= tape_pos)
	    {
	      __uint32 cur = *stack_pos;
	      stack_pos--;
 
	      if (datacenter[cur] != OPEN)
                  continue;
 
	      visited_count++;
 
	      datacenter[cur] = DFS_VISIT;
	      dfs_pos++;
	      *dfs_pos = cur;
 
	      __uint32 above = cur - a_width;
	      __uint32 left = cur - 1;
	      __uint32 right = cur + 1;
	      __uint32 below = cur + a_width;
 
	      __uint8 value, sum = 0;
 
	      stack_pos++;
 
              // A trick is being played here with 'value &= 1' so that
              // the following is straight-line code. (Conditional branches
              // will cause hiccups in the processor.) Since only open
              // rooms have value 1, they alone will be stored on the stack.
              // This avoids visiting closed or visited rooms.
 
	      value = datacenter[above];
              sum += value;
              value &= 1;
	      *stack_pos = above;
	      stack_pos += value;
 
	      value = datacenter[right];
              sum += value;
              value &= 1;
	      *stack_pos = right;
	      stack_pos += value;
 
	      value = datacenter[below];
              sum += value;
              value &= 1;
	      *stack_pos = below;
	      stack_pos += value;
 
	      value = datacenter[left];
              sum += value;
              value &= 1;
	      *stack_pos = left;
	      stack_pos += value;
 
              if (sum == DFS_VISIT && cur != end)
                 {
                   // If the sum of the values of the surrounding rooms
                   // is DFS_VISIT then it must be that all but one are
                   // closed - e.g. with p denoting the current position:
                   // 0 0 0
                   // 0 p 8
                   // 0 0 0
                   // This is a dead-end. Unless this position is also
                   // end, there will be no path through here that connects
                   // start and end.
                   // Set visited_count to max_steps only so the
                   // reachability check will fail.
 
                   visited_count = max_steps;
                   break;
                 }
 
	      stack_pos--;
	    }
 
          // Re-open those rooms which were marked DFS_VISIT within the
          // above loop.
 
	  for(; dfs_pos > dfs; dfs_pos--)
	    {
	      datacenter[*dfs_pos] = OPEN;
	    }
 
	  *tape_pos = (pos | VISIT_MASK);
	  datacenter[pos] = VISITED;
 
          // This is the reachability check. We will only append new
          // positions to visit on the tape if
 
	  if ((visited_count + steps) == max_steps)
	    {
	      __uint32 above = pos - a_width;
	      __uint32 left = pos - 1;
	      __uint32 right = pos + 1;
	      __uint32 below = pos + a_width;
 
	      __uint8 value;
 
	      tape_pos++;
 
              // This is using the same trick as the DFS loop for
              // reachability. There's no need to 'value &= 1'
              // because every room has either a 1 or 0 in it. If
              // we add a new position to the tape it is initially
              // added un-visited. We will mark it visited later
              // with VISIT_MASK.
 
	      value = datacenter[above];
	      *tape_pos = above;
	      tape_pos += value;
 
	      value = datacenter[right];
	      *tape_pos = right;
	      tape_pos += value;
 
	      value = datacenter[below];
	      *tape_pos = below;
	      tape_pos += value;
 
	      value = datacenter[left];
	      *tape_pos = left;
	      tape_pos += value;
 
	      tape_pos--;
	    }
 
          steps++;
	}
 
      // This is where the tape is used to rewind rooms occupied in the
      // datacenter. Above, on either sides of the if/else, was:
      //    *tape_pos = (pos | VISIT_MASK);
      // which set the high-bit in the tape position. Now, _bittestandreset,
      // returns that bit and sets it to 0. So the pseudocode reads:
      //    while the high bit at tape_pos is set
      //       clear the high bit
      //       mark the room at tape_pos as OPEN
      //       decrement tape_pos and steps
 
      while(_bittestandreset((long *)tape_pos, 31))
	{
	  datacenter[*tape_pos] = OPEN;
	  tape_pos--;
	  steps--;
	}
    }
 
  cout<<routes<<endl;
 
#ifdef MEASURETIME
  LARGE_INTEGER stop_time;
  QueryPerformanceCounter(&stop_time);
  LARGE_INTEGER frequency_time;
  QueryPerformanceFrequency(&frequency_time);
  double seconds = (double)(stop_time.QuadPart - start_time.QuadPart) /
    (double)(frequency_time.QuadPart);
  cout<<"Computed in "<<seconds<<" seconds."<<endl;
#endif MEASURETIME
}