[FrontPage] [TitleIndex] [WordIndex

Note: You are looking at a static copy of the former PineWiki site, used for class notes by James Aspnes from 2003 to 2012. Many mathematical formulas are broken, and there are likely to be other bugs as well. These will most likely not be fixed. You may be able to find more up-to-date versions of some of these notes at http://www.cs.yale.edu/homes/aspnes/#classes.

1. Searching a huge grid

For this assignment, you are to implement the A* search algorithm on a two-dimensional grid. Each position in the grid is represented by two coordinates of type int, packed into a struct position defined in search.h:

   1 struct position {
   2     int x;
   3     int y;
   4 };
   5 
   6 #define NO_PATH (-1)
   7 
   8 /* return length of shortest path from source to target 
   9  * using only points for which blocked returns 0,
  10  * or NO_PATH if there is no path.
  11  */
  12 int
  13 search(struct position source, struct position target, int (*blocked)(struct position));
search.h

The search routine takes a source and target position, together with a function blocked that returns nonzero for positions that are inaccessible. It should return the length of the shortest path from source to target that does not include any nodes for which blocked returns a nonzero value, where each new node in the path is directly north, east, south, or west of its predecessor. It does not need to compute the actual path.

If you can determine that there is no path (typically this will be because some region around source has only blocked neighbors), you should return NO_PATH, defined in search.h.

2. Examples

Here is a path of length 4 (marked with *) with no obstacles:

..t
..*
s**

Here is a path of length 8 that routes around some obstacles:

.****
**#.*
s.#.t
..#..
..#..
..#..

Note that neither of these paths is unique.

3. A* search

Because the grid is quite large (264 positions), you should probably not attempt to reuse your code from ../HW3. We recommend using A* search, a variant of Dijkstra's algorithm that uses a heuristic to favor expanding unexplored nodes that are closer to the target.

Recall that Dijkstra's algorithm keeps track of two sets: the closed nodes that have already been expanded and the open nodes that have been discovered but not expanded yet. Initially, the open nodes consist only of the source node. At each step, Dijkstra's algorithm chooses the open node with the smallest computed distance, marks it closed, and inserts all of its neighbors that are not already marked closed into the open set, possibly updating their distances if they are in the open set already. It continues until the target appears in the closed set. Dijkstra's algorithm has the property that, as long as there are no negative-weight edges, it has always computed the correct minimum distance to each node when that node leaves the open set.

For A*, we do the same thing, but instead of sorting nodes in the open set by d(source, node), we sort by d(source, node) + h(node, target), where h(node, target) is a lower bound on the distance from x to target. For the Taxicab geometry assumed in the problem, a good choice for h is "Manhattan distance", defined as |node.x - target.x| + |node.y - target.y|, the sum of the differences in the coordinates. This would be the length of a path uninterrupted by obstacles. The A* algorithm guarantees that it finds the shortest path as long as h is a non-negative lower bound on the actual distance, and that for any edge, the change in h is less than or equal to the length of the edge (which is true for this particular h and all edges in the grid). The reason for this is that under these conditions, we can pretend that each edge uv really has length l(uv)-h(u,target)+h(v,target), and then A* turns into Dijkstra's algorithm.

The advantage of A* over Dijkstra's algorithm is that it doesn't search through all nodes with d(source, node) < d(source, target). In the grid, this can produce quadratic speedups in the best case. There are many descriptions of A* on-line, including this visualization of the algorithm in action, a detailed WikiPedia entry, and these very detailed notes on implementing A* on a grid, aimed at game programmers.

4. What to submit

In addition to the search.h header file, we have provided a Makefile and test code testSearch.c. You should provide the missing search.c file that implements search. If you would like to structure your solution as multiple files, you should feel free to do so; in this case you should submit these additional files together with your own Makefile that will compile testSearch from your files and testSearch.c (which may be different in the final test script from the one we provide here).

Submit your files as usual with /c/cs223/bin/submit 10. You can run the public test script /c/cs223/Hwk10/test.public on your submitted files with /c/cs223/bin/testit 10.

5. Sample solution

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <assert.h>
   4 
   5 #include "search.h"
   6 
   7 /* we need a hash table to mark positions we have visited already */
   8 
   9 #define INITIAL_HASH_SIZE (1024)
  10 struct hash {
  11     int size;      /* number of slots in table */
  12     int count;     /* number of slots actually used */
  13     struct hashElt {
  14         int used;
  15         struct position pos;
  16     } *table;
  17 };
  18 
  19 static struct hash *
  20 hashCreate(int size)
  21 {
  22     struct hash *h;
  23     int i;
  24 
  25     h = malloc(sizeof(struct hash));
  26     h->size = size;
  27     h->count = 0;
  28     h->table = malloc(sizeof(struct hashElt) * h->size);
  29 
  30     for(i = 0; i < h->size; i++) {
  31         h->table[i].used = 0;
  32     }
  33 
  34     return h;
  35 }
  36 
  37 #define X_MULTIPLIER (19299141)
  38 #define Y_MULTIPLIER (97)
  39 
  40 static unsigned int
  41 hashPos(struct position pos)
  42 {
  43     return pos.x * X_MULTIPLIER + pos.y * Y_MULTIPLIER;
  44 }
  45 
  46 static void
  47 hashInsert(struct hash *h, struct position pos)
  48 {
  49     struct hash *h2;
  50     int i;
  51 
  52     if(h->count > h->size / 2) {
  53         /* double the size */
  54         h2 = hashCreate(h->size * 2);
  55         
  56         for(i = 0; i < h->size; i++) {
  57             if(h->table[i].used) {
  58                 hashInsert(h2, h->table[i].pos);
  59             }
  60         }
  61 
  62         free(h->table);
  63         *h = *h2;
  64 
  65         free(h2);
  66     }
  67 
  68     for(i = hashPos(pos) % h->size; h->table[i].used; i = (i+1) % h->size);
  69 
  70     h->table[i].used = 1;
  71     h->table[i].pos = pos;
  72     h->count++;
  73 }
  74 
  75 static int
  76 hashContains(struct hash *h, struct position pos)
  77 {
  78     int i;
  79 
  80     for(i = hashPos(pos) % h->size; h->table[i].used; i = (i+1) % h->size) {
  81         if(h->table[i].pos.x == pos.x && h->table[i].pos.y == pos.y) {
  82             return 1;
  83         }
  84     }
  85 
  86     return 0;
  87 }
  88 
  89 static void
  90 hashDestroy(struct hash *h)
  91 {
  92     free(h->table);
  93     free(h);
  94 }
  95 
  96 /* we need a heap to implement the A* priority queue */
  97 /* this prefers lower estimate */
  98 /* then breaks ties by greater distance */
  99 
 100 #define INITIAL_HEAP_SIZE (1024)
 101 
 102 #define Parent(x) ((x-1)/2)
 103 #define Child(x,i) (2*(x)+1+(i))
 104 
 105 struct heap {
 106     int size;   /* allocated positions in heap */
 107     int top;    /* first unused position in heap */
 108     struct heapElt {
 109         int estimate;   /* d(source, pos) + minimum possible d(pos, target) */
 110         int distance;   /* d(source, pos) */
 111         struct position pos;
 112     } *contents;
 113 };
 114 
 115 /* create an empty heap */
 116 static struct heap *
 117 heapCreate(void)
 118 {
 119     struct heap *h;
 120 
 121     h = malloc(sizeof(*h));
 122     assert(h);
 123 
 124     h->size = INITIAL_HEAP_SIZE;
 125     h->top = 0;
 126     h->contents = malloc(sizeof(struct heapElt) * h->size);
 127     assert(h->contents);
 128 
 129     return h;
 130 }
 131 
 132 static int
 133 heapLessThan(const struct heapElt *e1, const struct heapElt *e2)
 134 {
 135     return e1->estimate < e2->estimate;
 136 }
 137 
 138 static void
 139 heapSwap(struct heapElt *e1, struct heapElt *e2)
 140 {
 141     struct heapElt temp;
 142 
 143     temp = *e1;
 144     *e1 = *e2;
 145     *e2 = temp;
 146 }
 147 
 148 static void
 149 heapInsert(struct heap *h, struct heapElt new)
 150 {
 151     int i;
 152 
 153     if(h->top >= h->size) {
 154         h->size *= 2;
 155         h->contents = realloc(h->contents, sizeof(struct heapElt) * h->size);
 156         assert(h->contents);
 157     }
 158 
 159     i = h->top++;
 160 
 161     h->contents[i] = new;
 162 
 163     /* float up */
 164     while(i > 0 && heapLessThan(&h->contents[i], &h->contents[Parent(i)])) {
 165         heapSwap(&h->contents[i], &h->contents[Parent(i)]);
 166         i = Parent(i);
 167     }
 168 }
 169 
 170 static int
 171 heapIsEmpty(struct heap *h)
 172 {
 173     return h->top <= 0;
 174 }
 175 
 176 static struct heapElt
 177 heapDeleteMin(struct heap *h)
 178 {
 179     struct heapElt ret;
 180     int i;
 181     int smallerKid;
 182 
 183     assert(h->top > 0);
 184 
 185     ret = h->contents[0];
 186 
 187     h->contents[0] = h->contents[--(h->top)];
 188 
 189     i = 0;
 190 
 191     /* float down */
 192     while(Child(i,0) < h->top) {
 193         /* still have at least one kid */
 194         /* compute smaller kid */
 195         if(Child(i,1) >= h->top || heapLessThan(&h->contents[Child(i,0)], &h->contents[Child(i,1)])) {
 196             smallerKid = Child(i,0);
 197         } else {
 198             smallerKid = Child(i,1);
 199         }
 200 
 201         /* swap with smaller kid if kid is smaller */
 202         if(heapLessThan(&h->contents[smallerKid], &h->contents[i])) {
 203             heapSwap(&h->contents[smallerKid], &h->contents[i]);
 204             i = smallerKid;
 205         } else {
 206             break;
 207         }
 208     }
 209 
 210     return ret;
 211 }
 212 
 213 static void
 214 heapDestroy(struct heap *h)
 215 {
 216     free(h->contents);
 217     free(h);
 218 }
 219 
 220 /* actual search routine */
 221 
 222 static int
 223 iabs(int x)
 224 {
 225     return x < 0 ? -x : x;
 226 }
 227 
 228 static int
 229 distLowerBound(struct position source, struct position target)
 230 {
 231     return iabs(source.x - target.x) + iabs(source.y - target.y);
 232 }
 233 
 234 #define NUM_DIRS (4)
 235 
 236 static struct position Directions[NUM_DIRS] = {
 237     { -1, 0 },
 238     { 0, 1 },
 239     { 1, 0 },
 240     { 0, -1 }
 241 };
 242 
 243 int
 244 search(struct position source, struct position target, int (*blocked)(struct position))
 245 {
 246     struct hash *seen;    /* previously visited positions */
 247     struct heap *pq;      /* nodes to process */
 248     struct heapElt elt; 
 249     struct heapElt neighbor; 
 250     int dir;
 251 
 252     if(blocked(source) || blocked(target)) {
 253         return NO_PATH;
 254     }
 255 
 256     seen = hashCreate(INITIAL_HASH_SIZE);
 257     pq = heapCreate();
 258 
 259     /* push source to start with */
 260     elt.pos = source;
 261     elt.distance = 0;
 262     elt.estimate = distLowerBound(source, target);
 263 
 264     heapInsert(pq, elt);
 265 
 266     while(!heapIsEmpty(pq)) {
 267         elt = heapDeleteMin(pq);
 268 
 269         if(hashContains(seen, elt.pos)) continue; /* saw a closer copy */
 270 
 271         /* else */
 272         hashInsert(seen, elt.pos);
 273 
 274         if(elt.pos.x == target.x && elt.pos.y == target.y) {
 275             /* we are done */
 276             hashDestroy(seen);
 277             heapDestroy(pq);
 278 
 279             return elt.distance;
 280         }
 281 
 282         for(dir = 0; dir < NUM_DIRS; dir++) {
 283             neighbor.pos.x = elt.pos.x + Directions[dir].x;
 284             neighbor.pos.y = elt.pos.y + Directions[dir].y;
 285 
 286             if(!blocked(neighbor.pos) && !hashContains(seen, neighbor.pos)) {
 287                 /* push it */
 288                 neighbor.distance = elt.distance + 1;
 289                 neighbor.estimate = neighbor.distance + distLowerBound(neighbor.pos, target);
 290 
 291                 heapInsert(pq, neighbor);
 292             }
 293         }
 294     }
 295 
 296     heapDestroy(pq);
 297     hashDestroy(seen);
 298 
 299     return NO_PATH;   /* no path */
 300 }
search.c

2014-06-17 11:58