src/autoplace.c

   1 /* $Id$ */
   2
   3 /*
   4  *                            COPYRIGHT
   5  *
   6  *  PCB, interactive printed circuit board design
   7  *  Copyright (C) 1994,1995,1996 Thomas Nau
   8  *  Copyright (C) 1998,1999,2000,2001 harry eaton
   9  *
  10  *  This program is free software; you can redistribute it and/or modify
  11  *  it under the terms of the GNU General Public License as published by
  12  *  the Free Software Foundation; either version 2 of the License, or
  13  *  (at your option) any later version.
  14  *
  15  *  This program is distributed in the hope that it will be useful,
  16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  *  GNU General Public License for more details.
  19  *
  20  *  You should have received a copy of the GNU General Public License
  21  *  along with this program; if not, write to the Free Software
  22  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  23  *
  24  *  Contact addresses for paper mail and Email:
  25  *  harry eaton, 6697 Buttonhole Ct, Columbia, MD 21044 USA
  26  *  haceaton@aplcomm.jhuapl.edu
  27  *
  28  */
  29
  30 /*
  31  * This moduel, autoplace.c, was written by and is
  32  * Copyright (c) 2001 C. Scott Ananian
  33  */
  34
  35 /* functions used to autoplace elements.
  36  */
  37
  38 #ifdef HAVE_CONFIG_H
  39 #include "config.h"
  40 #endif
  41
  42 #include <assert.h>
  43 #include <math.h>
  44 #include <memory.h>
  45 #include <stdlib.h>
  46
  47 #include "global.h"
  48
  49 #include "autoplace.h"
  50 #include "box.h"
  51 #include "compat.h"
  52 #include "data.h"
  53 #include "draw.h"
  54 #include "error.h"
  55 #include "intersect.h"
  56 #include "rtree.h"
  57 #include "macro.h"
  58 #include "mirror.h"
  59 #include "misc.h"
  60 #include "move.h"
  61 #include "mymem.h"
  62 #include "rats.h"
  63 #include "remove.h"
  64 #include "rotate.h"
  65
  66 #ifdef HAVE_LIBDMALLOC
  67 #include <dmalloc.h>
  68 #endif
  69
  70 RCSID ("$Id$");
  71
  72 #define EXPANDRECTXY(r1, x1, y1, x2, y2) { \
  73   r1->X1=MIN(r1->X1, x1); r1->Y1=MIN(r1->Y1, y1); \
  74   r1->X2=MAX(r1->X2, x2); r1->Y2=MAX(r1->Y2, y2); \
  75 }
  76 #define EXPANDRECT(r1, r2) EXPANDRECTXY(r1, r2->X1, r2->Y1, r2->X2, r2->Y2)
  77
  78 /* ---------------------------------------------------------------------------
  79  * some local prototypes
  80  */
  81 static double ComputeCost (NetListTypePtr Nets, double T0, double T);
  82
  83 /* ---------------------------------------------------------------------------
  84  * some local types
  85  */
  86 const struct
  87 {
  88   float via_cost;
  89   float congestion_penalty;     /* penalty length / unit area */
  90   float overlap_penalty_min;    /* penalty length / unit area at start */
  91   float overlap_penalty_max;    /* penalty length / unit area at end */
  92   float out_of_bounds_penalty;  /* assessed for each component oob */
  93   float overall_area_penalty;   /* penalty length / unit area */
  94   float matching_neighbor_bonus;        /* length bonus per same-type neigh. */
  95   float aligned_neighbor_bonus; /* length bonus per aligned neigh. */
  96   float oriented_neighbor_bonus;        /* length bonus per same-rot neigh. */
  97 #if 0
  98   float pin_alignment_bonus;    /* length bonus per exact alignment */
  99   float bound_alignment_bonus;  /* length bonus per exact alignment */
 100 #endif
 101   float m;                      /* annealing stage cutoff constant */
 102   float gamma;                  /* annealing schedule constant */
 103   int good_ratio;               /* ratio of moves to good moves for halting */
 104   bool fast;                    /* ignore SMD/pin conflicts */
 105   int large_grid_size;          /*snap perturbations to this grid when T is high */
 106   int small_grid_size;          /* snap to this grid when T is small. */
 107 }
 108 /* wire cost is manhattan distance (in mils), thus 1 inch = 1000 */
 109 CostParameter =
 110 {
 111   3e3,                          /* via cost */
 112     2e-2,                       /* congestion penalty */
 113     1e-2,                       /* initial overlap penalty */
 114     1e2,                        /* final overlap penalty */
 115     1e3,                        /* out of bounds penalty */
 116     1e0,                        /* penalty for total area used */
 117     1e0,                        /* subtract 1000 from cost for every same-type neighbor */
 118     1e0,                        /* subtract 1000 from cost for every aligned neighbor */
 119     1e0,                        /* subtract 1000 from cost for every same-rotation neighbor */
 120     20,                         /* move on when each module has been profitably moved 20 times */
 121     0.75,                       /* annealing schedule constant: 0.85 */
 122     40,                         /* halt when there are 60 times as many moves as good moves */
 123     false,                      /* don't ignore SMD/pin conflicts */
 124     100,                        /* coarse grid is 100 mils */
 125     10,                         /* fine grid is 10 mils */
 126 };
 127
 128 typedef struct
 129 {
 130   ElementTypePtr *element;
 131   Cardinal elementN;
 132 }
 133 ElementPtrListType;
 134
 135 typedef struct
 136 {
 137   ElementTypePtr element;
 138   enum
 139   { SHIFT, ROTATE, EXCHANGE }
 140   which;
 141   LocationType DX, DY;          /* for shift */
 142   BYTE rotate;                  /* for rotate/flip */
 143   ElementTypePtr other;         /* for exchange */
 144 }
 145 PerturbationType;
 146
 147 /* ---------------------------------------------------------------------------
 148  * some local identifiers
 149  */
 150
 151 /* ---------------------------------------------------------------------------
 152  * Update the X, Y and group position information stored in the NetList after
 153  * elements have possibly been moved, rotated, flipped, etc.
 154  */
 155 static void
 156 UpdateXY (NetListTypePtr Nets)
 157 {
 158   Cardinal SLayer, CLayer;
 159   Cardinal i, j;
 160   /* find layer groups of the component side and solder side */
 161   SLayer = GetLayerGroupNumberByNumber (max_layer + SOLDER_LAYER);
 162   CLayer = GetLayerGroupNumberByNumber (max_layer + COMPONENT_LAYER);
 163   /* update all nets */
 164   for (i = 0; i < Nets->NetN; i++)
 165     {
 166       for (j = 0; j < Nets->Net[i].ConnectionN; j++)
 167         {
 168           ConnectionTypePtr c = &(Nets->Net[i].Connection[j]);
 169           switch (c->type)
 170             {
 171             case PAD_TYPE:
 172               c->group = TEST_FLAG (ONSOLDERFLAG,
 173                                     (ElementTypePtr) c->ptr1)
 174                 ? SLayer : CLayer;
 175               c->X = ((PadTypePtr) c->ptr2)->Point1.X;
 176               c->Y = ((PadTypePtr) c->ptr2)->Point1.Y;
 177               break;
 178             case PIN_TYPE:
 179               c->group = SLayer;        /* any layer will do */
 180               c->X = ((PinTypePtr) c->ptr2)->X;
 181               c->Y = ((PinTypePtr) c->ptr2)->Y;
 182               break;
 183             default:
 184               Message ("Odd connection type encountered in " "UpdateXY");
 185               break;
 186             }
 187         }
 188     }
 189 }
 190
 191 /* ---------------------------------------------------------------------------
 192  * Create a list of selected elements.
 193  */
 194 static PointerListType
 195 collectSelectedElements ()
 196 {
 197   PointerListType list = { 0, 0, NULL };
 198   ELEMENT_LOOP (PCB->Data);
 199   {
 200     if (TEST_FLAG (SELECTEDFLAG, element))
 201       {
 202         ElementTypePtr *epp = (ElementTypePtr *) GetPointerMemory (&list);
 203         *epp = element;
 204       }
 205   }
 206   END_LOOP;
 207   return list;
 208 }
 209
 210 #if 0                           /* only for debugging box lists */
 211 #include "create.h"
 212 /* makes a line on the solder layer surrounding all boxes in blist */
 213 static void
 214 showboxes (BoxListTypePtr blist)
 215 {
 216   Cardinal i;
 217   LayerTypePtr SLayer = &(PCB->Data->Layer[max_layer + SOLDER_LAYER]);
 218   for (i = 0; i < blist->BoxN; i++)
 219     {
 220       CreateNewLineOnLayer (SLayer, blist->Box[i].X1, blist->Box[i].Y1,
 221                             blist->Box[i].X2, blist->Box[i].Y1, 1, 1, 0);
 222       CreateNewLineOnLayer (SLayer, blist->Box[i].X1, blist->Box[i].Y2,
 223                             blist->Box[i].X2, blist->Box[i].Y2, 1, 1, 0);
 224       CreateNewLineOnLayer (SLayer, blist->Box[i].X1, blist->Box[i].Y1,
 225                             blist->Box[i].X1, blist->Box[i].Y2, 1, 1, 0);
 226       CreateNewLineOnLayer (SLayer, blist->Box[i].X2, blist->Box[i].Y1,
 227                             blist->Box[i].X2, blist->Box[i].Y2, 1, 1, 0);
 228     }
 229 }
 230 #endif
 231
 232 /* ---------------------------------------------------------------------------
 233  * Helper function to compute "closest neighbor" for a box in a rtree.
 234  * The closest neighbor on a certain side is the closest one in a trapezoid
 235  * emanating from that side.
 236  */
 237 /*------ r_find_neighbor ------*/
 238 struct r_neighbor_info
 239 {
 240   const BoxType *neighbor;
 241   BoxType trap;
 242   direction_t search_dir;
 243 };
 244 #define ROTATEBOX(box) { LocationType t;\
 245     t = (box).X1; (box).X1 = - (box).Y1; (box).Y1 = t;\
 246     t = (box).X2; (box).X2 = - (box).Y2; (box).Y2 = t;\
 247     t = (box).X1; (box).X1 =   (box).X2; (box).X2 = t;\
 248 }
 249 /* helper methods for __r_find_neighbor */
 250 static int
 251 __r_find_neighbor_reg_in_sea (const BoxType * region, void *cl)
 252 {
 253   struct r_neighbor_info *ni = (struct r_neighbor_info *) cl;
 254   BoxType query = *region;
 255   ROTATEBOX_TO_NORTH (query, ni->search_dir);
 256   /*  ______________ __ trap.y1     __
 257    *  \            /               |__| query rect.
 258    *   \__________/  __ trap.y2
 259    *   |          |
 260    *   trap.x1    trap.x2   sides at 45-degree angle
 261    */
 262   return (query.Y2 > ni->trap.Y1) && (query.Y1 < ni->trap.Y2) &&
 263     (query.X2 + ni->trap.Y2 > ni->trap.X1 + query.Y1) &&
 264     (query.X1 + query.Y1 < ni->trap.X2 + ni->trap.Y2);
 265 }
 266 static int
 267 __r_find_neighbor_rect_in_reg (const BoxType * box, void *cl)
 268 {
 269   struct r_neighbor_info *ni = (struct r_neighbor_info *) cl;
 270   BoxType query = *box;
 271   int r;
 272   ROTATEBOX_TO_NORTH (query, ni->search_dir);
 273   /*  ______________ __ trap.y1     __
 274    *  \            /               |__| query rect.
 275    *   \__________/  __ trap.y2
 276    *   |          |
 277    *   trap.x1    trap.x2   sides at 45-degree angle
 278    */
 279   r = (query.Y2 > ni->trap.Y1) && (query.Y1 < ni->trap.Y2) &&
 280     (query.X2 + ni->trap.Y2 > ni->trap.X1 + query.Y1) &&
 281     (query.X1 + query.Y1 < ni->trap.X2 + ni->trap.Y2);
 282   r = r && (query.Y2 <= ni->trap.Y2);
 283   if (r)
 284     {
 285       ni->trap.Y1 = query.Y2;
 286       ni->neighbor = box;
 287     }
 288   return r;
 289 }
 290
 291 /* main r_find_neighbor routine.  Returns NULL if no neighbor in the
 292  * requested direction. */
 293 static const BoxType *
 294 r_find_neighbor (rtree_t * rtree, const BoxType * box,
 295                  direction_t search_direction)
 296 {
 297   struct r_neighbor_info ni;
 298   BoxType bbox;
 299
 300   ni.neighbor = NULL;
 301   ni.trap = *box;
 302   ni.search_dir = search_direction;
 303
 304   bbox.X1 = bbox.Y1 = 0;
 305   bbox.X2 = PCB->MaxWidth;
 306   bbox.Y2 = PCB->MaxHeight;
 307   /* rotate so that we can use the 'north' case for everything */
 308   ROTATEBOX_TO_NORTH (bbox, search_direction);
 309   ROTATEBOX_TO_NORTH (ni.trap, search_direction);
 310   /* shift Y's such that trap contains full bounds of trapezoid */
 311   ni.trap.Y2 = ni.trap.Y1;
 312   ni.trap.Y1 = bbox.Y1;
 313   /* do the search! */
 314   r_search (rtree, NULL,
 315             __r_find_neighbor_reg_in_sea, __r_find_neighbor_rect_in_reg, &ni);
 316   return ni.neighbor;
 317 }
 318
 319 /* ---------------------------------------------------------------------------
 320  * Compute cost function.
 321  *  note that area overlap cost is correct for SMD devices: SMD devices on
 322  *  opposite sides of the board don't overlap.
 323  *
 324  * Algorithms follow those described in sections 4.1 of
 325  *  "Placement and Routing of Electronic Modules" edited by Michael Pecht
 326  *  Marcel Dekker, Inc. 1993.  ISBN: 0-8247-8916-4 TK7868.P7.P57 1993
 327  */
 328 static double
 329 ComputeCost (NetListTypePtr Nets, double T0, double T)
 330 {
 331   double W = 0;                 /* wire cost */
 332   double delta1 = 0;            /* wire congestion penalty function */
 333   double delta2 = 0;            /* module overlap penalty function */
 334   double delta3 = 0;            /* out of bounds penalty */
 335   double delta4 = 0;            /* alignment bonus */
 336   double delta5 = 0;            /* total area penalty */
 337   Cardinal i, j;
 338   LocationType minx, maxx, miny, maxy;
 339   bool allpads, allsameside;
 340   Cardinal thegroup;
 341   BoxListType bounds = { 0, 0, NULL };  /* save bounding rectangles here */
 342   BoxListType solderside = { 0, 0, NULL };      /* solder side component bounds */
 343   BoxListType componentside = { 0, 0, NULL };   /* component side bounds */
 344   /* make sure the NetList have the proper updated X and Y coords */
 345   UpdateXY (Nets);
 346   /* wire length term.  approximated by half-perimeter of minimum
 347    * rectangle enclosing the net.  Note that we penalize vias in
 348    * all-SMD nets by making the rectangle a cube and weighting
 349    * the "layer height" of the net. */
 350   for (i = 0; i < Nets->NetN; i++)
 351     {
 352       NetTypePtr n = &Nets->Net[i];
 353       if (n->ConnectionN < 2)
 354         continue;               /* no cost to go nowhere */
 355       minx = maxx = n->Connection[0].X;
 356       miny = maxy = n->Connection[0].Y;
 357       thegroup = n->Connection[0].group;
 358       allpads = (n->Connection[0].type == PAD_TYPE);
 359       allsameside = true;
 360       for (j = 1; j < n->ConnectionN; j++)
 361         {
 362           ConnectionTypePtr c = &(n->Connection[j]);
 363           MAKEMIN (minx, c->X);
 364           MAKEMAX (maxx, c->X);
 365           MAKEMIN (miny, c->Y);
 366           MAKEMAX (maxy, c->Y);
 367           if (c->type != PAD_TYPE)
 368             allpads = false;
 369           if (c->group != thegroup)
 370             allsameside = false;
 371         }
 372       /* save bounding rectangle */
 373       {
 374         BoxTypePtr box = GetBoxMemory (&bounds);
 375         box->X1 = minx;
 376         box->Y1 = miny;
 377         box->X2 = maxx;
 378         box->Y2 = maxy;
 379       }
 380       /* okay, add half-perimeter to cost! */
 381       W += (maxx - minx) / 100 + (maxy - miny) / 100 +
 382         ((allpads && !allsameside) ? CostParameter.via_cost : 0);
 383     }
 384   /* now compute penalty function Wc which is proportional to
 385    * amount of overlap and congestion. */
 386   /* delta1 is congestion penalty function */
 387   delta1 = CostParameter.congestion_penalty *
 388     sqrt (fabs (ComputeIntersectionArea (&bounds)));
 389 #if 0
 390   printf ("Wire Congestion Area: %f\n", ComputeIntersectionArea (&bounds));
 391 #endif
 392   /* free bounding rectangles */
 393   FreeBoxListMemory (&bounds);
 394   /* now collect module areas (bounding rect of pins/pads) */
 395   /* two lists for solder side / component side. */
 396
 397   ELEMENT_LOOP (PCB->Data);
 398   {
 399     BoxListTypePtr thisside;
 400     BoxListTypePtr otherside;
 401     BoxTypePtr box;
 402     BoxTypePtr lastbox = NULL;
 403     BDimension thickness;
 404     BDimension clearance;
 405     if (TEST_FLAG (ONSOLDERFLAG, element))
 406       {
 407         thisside = &solderside;
 408         otherside = &componentside;
 409       }
 410     else
 411       {
 412         thisside = &componentside;
 413         otherside = &solderside;
 414       }
 415     box = GetBoxMemory (thisside);
 416     /* protect against elements with no pins/pads */
 417     if (element->PinN == 0 && element->PadN == 0)
 418       continue;
 419     /* initialize box so that it will take the dimensions of
 420      * the first pin/pad */
 421     box->X1 = MAX_COORD;
 422     box->Y1 = MAX_COORD;
 423     box->X2 = -MAX_COORD;
 424     box->Y2 = -MAX_COORD;
 425     PIN_LOOP (element);
 426     {
 427       thickness = pin->Thickness / 2;
 428       clearance = pin->Clearance * 2;
 429     EXPANDRECTXY (box,
 430                     pin->X - (thickness + clearance),
 431                     pin->Y - (thickness + clearance),
 432                     pin->X + (thickness + clearance),
 433                     pin->Y + (thickness + clearance))}
 434     END_LOOP;
 435     PAD_LOOP (element);
 436     {
 437       thickness = pad->Thickness / 2;
 438       clearance = pad->Clearance * 2;
 439     EXPANDRECTXY (box,
 440                     MIN (pad->Point1.X,
 441                            pad->Point2.X) - (thickness +
 442                                                clearance),
 443                     MIN (pad->Point1.Y,
 444                            pad->Point2.Y) - (thickness +
 445                                                clearance),
 446                     MAX (pad->Point1.X,
 447                            pad->Point2.X) + (thickness +
 448                                                clearance),
 449                     MAX (pad->Point1.Y,
 450                            pad->Point2.Y) + (thickness + clearance))}
 451     END_LOOP;
 452     /* add a box for each pin to the "opposite side":
 453      * surface mount components can't sit on top of pins */
 454     if (!CostParameter.fast)
 455       PIN_LOOP (element);
 456     {
 457       box = GetBoxMemory (otherside);
 458       thickness = pin->Thickness / 2;
 459       clearance = pin->Clearance * 2;
 460       /* we ignore clearance here */
 461       /* (otherwise pins don't fit next to each other) */
 462       box->X1 = pin->X - thickness;
 463       box->Y1 = pin->Y - thickness;
 464       box->X2 = pin->X + thickness;
 465       box->Y2 = pin->Y + thickness;
 466       /* speed hack! coalesce with last box if we can */
 467       if (lastbox != NULL &&
 468           ((lastbox->X1 == box->X1 &&
 469             lastbox->X2 == box->X2 &&
 470             MIN (abs (lastbox->Y1 - box->Y2),
 471                  abs (box->Y1 - lastbox->Y2)) <
 472             clearance) || (lastbox->Y1 == box->Y1
 473                            && lastbox->Y2 == box->Y2
 474                            &&
 475                            MIN (abs
 476                                 (lastbox->X1 -
 477                                  box->X2),
 478                                 abs (box->X1 - lastbox->X2)) < clearance)))
 479         {
 480           EXPANDRECT (lastbox, box);
 481           otherside->BoxN--;
 482         }
 483       else
 484         lastbox = box;
 485     }
 486     END_LOOP;
 487     /* assess out of bounds penalty */
 488     if (element->VBox.X1 < 0 ||
 489         element->VBox.Y1 < 0 ||
 490         element->VBox.X2 > PCB->MaxWidth || element->VBox.Y2 > PCB->MaxHeight)
 491       delta3 += CostParameter.out_of_bounds_penalty;
 492   }
 493   END_LOOP;
 494   /* compute intersection area of module areas box list */
 495   delta2 = sqrt (fabs (ComputeIntersectionArea (&solderside) +
 496                        ComputeIntersectionArea (&componentside))) *
 497     (CostParameter.overlap_penalty_min +
 498      (1 - (T / T0)) * CostParameter.overlap_penalty_max);
 499 #if 0
 500   printf ("Module Overlap Area (solder): %f\n",
 501           ComputeIntersectionArea (&solderside));
 502   printf ("Module Overlap Area (component): %f\n",
 503           ComputeIntersectionArea (&componentside));
 504 #endif
 505   FreeBoxListMemory (&solderside);
 506   FreeBoxListMemory (&componentside);
 507   /* reward pin/pad x/y alignment */
 508   /* score higher if pins/pads belong to same *type* of component */
 509   /* XXX: subkey should be *distance* from thing aligned with, so that
 510    * aligning to something far away isn't profitable */
 511   {
 512     /* create r tree */
 513     PointerListType seboxes = { 0, 0, NULL }
 514     , ceboxes =
 515     {
 516     0, 0, NULL};
 517     struct ebox
 518     {
 519       BoxType box;
 520       ElementTypePtr element;
 521     };
 522     direction_t dir[4] = { NORTH, EAST, SOUTH, WEST };
 523     struct ebox **boxpp, *boxp;
 524     rtree_t *rt_s, *rt_c;
 525     int factor;
 526     ELEMENT_LOOP (PCB->Data);
 527     {
 528       boxpp = (struct ebox **)
 529         GetPointerMemory (TEST_FLAG (ONSOLDERFLAG, element) ?
 530                           &seboxes : &ceboxes);
 531       *boxpp = malloc (sizeof (**boxpp));
 532       if (*boxpp == NULL )
 533         {
 534           fprintf (stderr, "malloc() failed in %s\n", __FUNCTION__);
 535           exit (1);
 536         }
 537
 538       (*boxpp)->box = element->VBox;
 539       (*boxpp)->element = element;
 540     }
 541     END_LOOP;
 542     rt_s = r_create_tree ((const BoxType **) seboxes.Ptr, seboxes.PtrN, 1);
 543     rt_c = r_create_tree ((const BoxType **) ceboxes.Ptr, ceboxes.PtrN, 1);
 544     FreePointerListMemory (&seboxes);
 545     FreePointerListMemory (&ceboxes);
 546     /* now, for each element, find its neighbor on all four sides */
 547     delta4 = 0;
 548     for (i = 0; i < 4; i++)
 549       ELEMENT_LOOP (PCB->Data);
 550     {
 551       boxp = (struct ebox *)
 552         r_find_neighbor (TEST_FLAG (ONSOLDERFLAG, element) ?
 553                          rt_s : rt_c, &element->VBox, dir[i]);
 554       /* score bounding box alignments */
 555       if (!boxp)
 556         continue;
 557       factor = 1;
 558       if (element->Name[0].TextString &&
 559           boxp->element->Name[0].TextString &&
 560           0 == NSTRCMP (element->Name[0].TextString,
 561                         boxp->element->Name[0].TextString))
 562         {
 563           delta4 += CostParameter.matching_neighbor_bonus;
 564           factor++;
 565         }
 566       if (element->Name[0].Direction == boxp->element->Name[0].Direction)
 567         delta4 += factor * CostParameter.oriented_neighbor_bonus;
 568       if (element->VBox.X1 ==
 569           boxp->element->VBox.X1 ||
 570           element->VBox.X1 ==
 571           boxp->element->VBox.X2 ||
 572           element->VBox.X2 ==
 573           boxp->element->VBox.X1 ||
 574           element->VBox.X2 ==
 575           boxp->element->VBox.X2 ||
 576           element->VBox.Y1 ==
 577           boxp->element->VBox.Y1 ||
 578           element->VBox.Y1 ==
 579           boxp->element->VBox.Y2 ||
 580           element->VBox.Y2 ==
 581           boxp->element->VBox.Y1 ||
 582           element->VBox.Y2 == boxp->element->VBox.Y2)
 583         delta4 += factor * CostParameter.aligned_neighbor_bonus;
 584     }
 585     END_LOOP;
 586     /* free k-d tree memory */
 587     r_destroy_tree (&rt_s);
 588     r_destroy_tree (&rt_c);
 589   }
 590   /* penalize total area used by this layout */
 591   {
 592     LocationType minX = MAX_COORD, minY = MAX_COORD;
 593     LocationType maxX = -MAX_COORD, maxY = -MAX_COORD;
 594     ELEMENT_LOOP (PCB->Data);
 595     {
 596       MAKEMIN (minX, element->VBox.X1);
 597       MAKEMIN (minY, element->VBox.Y1);
 598       MAKEMAX (maxX, element->VBox.X2);
 599       MAKEMAX (maxY, element->VBox.Y2);
 600     }
 601     END_LOOP;
 602     if (minX < maxX && minY < maxY)
 603       delta5 = CostParameter.overall_area_penalty *
 604         sqrt ((double) (maxX - minX) * (maxY - minY) * 0.0001);
 605   }
 606   if (T == 5)
 607     {
 608       T = W + delta1 + delta2 + delta3 - delta4 + delta5;
 609       printf ("cost components are %.3f %.3f %.3f %.3f %.3f %.3f\n",
 610               W / T, delta1 / T, delta2 / T, delta3 / T, -delta4 / T,
 611               delta5 / T);
 612     }
 613   /* done! */
 614   return W + (delta1 + delta2 + delta3 - delta4 + delta5);
 615 }
 616
 617 /* ---------------------------------------------------------------------------
 618  * Perturb:
 619  *  1) flip SMD from solder side to component side or vice-versa.
 620  *  2) rotate component 90, 180, or 270 degrees.
 621  *  3) shift component random + or - amount in random direction.
 622  *     (magnitude of shift decreases over time)
 623  *  -- Only perturb selected elements (need count/list of selected?) --
 624  */
 625 PerturbationType
 626 createPerturbation (PointerListTypePtr selected, double T)
 627 {
 628   PerturbationType pt = { 0 };
 629   /* pick element to perturb */
 630   pt.element = (ElementTypePtr) selected->Ptr[random () % selected->PtrN];
 631   /* exchange, flip/rotate or shift? */
 632   switch (random () % ((selected->PtrN > 1) ? 3 : 2))
 633     {
 634     case 0:
 635       {                         /* shift! */
 636         int grid;
 637         double scaleX = MAX (250, MIN (sqrt (T), PCB->MaxWidth / 3));
 638         double scaleY = MAX (250, MIN (sqrt (T), PCB->MaxHeight / 3));
 639         pt.which = SHIFT;
 640         pt.DX = scaleX * 2 * ((((double) random ()) / RAND_MAX) - 0.5);
 641         pt.DY = scaleY * 2 * ((((double) random ()) / RAND_MAX) - 0.5);
 642         /* snap to grid. different grids for "high" and "low" T */
 643         grid = (T > 1000) ? CostParameter.large_grid_size :
 644           CostParameter.small_grid_size;
 645         /* (round away from zero) */
 646         pt.DX = ((pt.DX / grid) + SGN (pt.DX)) * grid;
 647         pt.DY = ((pt.DY / grid) + SGN (pt.DY)) * grid;
 648         /* limit DX/DY so we don't fall off board */
 649         pt.DX = MAX (pt.DX, -pt.element->VBox.X1);
 650         pt.DX = MIN (pt.DX, PCB->MaxWidth - pt.element->VBox.X2);
 651         pt.DY = MAX (pt.DY, -pt.element->VBox.Y1);
 652         pt.DY = MIN (pt.DY, PCB->MaxHeight - pt.element->VBox.Y2);
 653         /* all done but the movin' */
 654         break;
 655       }
 656     case 1:
 657       {                         /* flip/rotate! */
 658         /* only flip if it's an SMD component */
 659         bool isSMD = pt.element->PadN != 0;
 660         pt.which = ROTATE;
 661         pt.rotate = isSMD ? (random () & 3) : (1 + (random () % 3));
 662         /* 0 - flip; 1-3, rotate. */
 663         break;
 664       }
 665     case 2:
 666       {                         /* exchange! */
 667         pt.which = EXCHANGE;
 668         pt.other = (ElementTypePtr)
 669           selected->Ptr[random () % (selected->PtrN - 1)];
 670         if (pt.other == pt.element)
 671           pt.other = (ElementTypePtr) selected->Ptr[selected->PtrN - 1];
 672         /* don't allow exchanging a solderside-side SMD component
 673          * with a non-SMD component. */
 674         if ((pt.element->PinN != 0 /* non-SMD */  &&
 675              TEST_FLAG (ONSOLDERFLAG, pt.other)) ||
 676             (pt.other->PinN != 0 /* non-SMD */  &&
 677              TEST_FLAG (ONSOLDERFLAG, pt.element)))
 678           return createPerturbation (selected, T);
 679         break;
 680       }
 681     default:
 682       assert (0);
 683     }
 684   return pt;
 685 }
 686
 687 void
 688 doPerturb (PerturbationType * pt, bool undo)
 689 {
 690   LocationType bbcx, bbcy;
 691   /* compute center of element bounding box */
 692   bbcx = (pt->element->VBox.X1 + pt->element->VBox.X2) / 2;
 693   bbcy = (pt->element->VBox.Y1 + pt->element->VBox.Y2) / 2;
 694   /* do exchange, shift or flip/rotate */
 695   switch (pt->which)
 696     {
 697     case SHIFT:
 698       {
 699         LocationType DX = pt->DX, DY = pt->DY;
 700         if (undo)
 701           {
 702             DX = -DX;
 703             DY = -DY;
 704           }
 705         MoveElementLowLevel (PCB->Data, pt->element, DX, DY);
 706         return;
 707       }
 708     case ROTATE:
 709       {
 710         BYTE b = pt->rotate;
 711         if (undo)
 712           b = (4 - b) & 3;
 713         /* 0 - flip; 1-3, rotate. */
 714         if (b)
 715           RotateElementLowLevel (PCB->Data, pt->element, bbcx, bbcy, b);
 716         else
 717           {
 718             LocationType y = pt->element->VBox.Y1;
 719             MirrorElementCoordinates (PCB->Data, pt->element, 0);
 720             /* mirroring moves the element.  move it back. */
 721             MoveElementLowLevel (PCB->Data, pt->element, 0,
 722                                  y - pt->element->VBox.Y1);
 723           }
 724         return;
 725       }
 726     case EXCHANGE:
 727       {
 728         /* first exchange positions */
 729         LocationType x1 = pt->element->VBox.X1;
 730         LocationType y1 = pt->element->VBox.Y1;
 731         LocationType x2 = pt->other->BoundingBox.X1;
 732         LocationType y2 = pt->other->BoundingBox.Y1;
 733         MoveElementLowLevel (PCB->Data, pt->element, x2 - x1, y2 - y1);
 734         MoveElementLowLevel (PCB->Data, pt->other, x1 - x2, y1 - y2);
 735         /* then flip both elements if they are on opposite sides */
 736         if (TEST_FLAG (ONSOLDERFLAG, pt->element) !=
 737             TEST_FLAG (ONSOLDERFLAG, pt->other))
 738           {
 739             PerturbationType mypt;
 740             mypt.element = pt->element;
 741             mypt.which = ROTATE;
 742             mypt.rotate = 0;    /* flip */
 743             doPerturb (&mypt, undo);
 744             mypt.element = pt->other;
 745             doPerturb (&mypt, undo);
 746           }
 747         /* done */
 748         return;
 749       }
 750     default:
 751       assert (0);
 752     }
 753 }
 754
 755 /* ---------------------------------------------------------------------------
 756  * Auto-place selected components.
 757  */
 758 bool
 759 AutoPlaceSelected (void)
 760 {
 761   NetListTypePtr Nets;
 762   PointerListType Selected = { 0, 0, NULL };
 763   PerturbationType pt;
 764   double C0, T0;
 765   bool changed = false;
 766
 767   /* (initial netlist processing copied from AddAllRats) */
 768   /* the netlist library has the text form
 769    * ProcNetlist fills in the Netlist
 770    * structure the way the final routing
 771    * is supposed to look
 772    */
 773   Nets = ProcNetlist (&PCB->NetlistLib);
 774   if (!Nets)
 775     {
 776       Message (_("Can't add rat lines because no netlist is loaded.\n"));
 777       goto done;
 778     }
 779
 780   Selected = collectSelectedElements ();
 781   if (Selected.PtrN == 0)
 782     {
 783       Message (_("No elements selected to autoplace.\n"));
 784       goto done;
 785     }
 786
 787   /* simulated annealing */
 788   {                             /* compute T0 by doing a random series of moves. */
 789     const int TRIALS = 10;
 790     const double Tx = 3e5, P = 0.95;
 791     double Cs = 0.0;
 792     int i;
 793     C0 = ComputeCost (Nets, Tx, Tx);
 794     for (i = 0; i < TRIALS; i++)
 795       {
 796         pt = createPerturbation (&Selected, 1e6);
 797         doPerturb (&pt, false);
 798         Cs += fabs (ComputeCost (Nets, Tx, Tx) - C0);
 799         doPerturb (&pt, true);
 800       }
 801     T0 = -(Cs / TRIALS) / log (P);
 802     printf ("Initial T: %f\n", T0);
 803   }
 804   /* now anneal in earnest */
 805   {
 806     double T = T0;
 807     long steps = 0;
 808     int good_moves = 0, moves = 0;
 809     const int good_move_cutoff = CostParameter.m * Selected.PtrN;
 810     const int move_cutoff = 2 * good_move_cutoff;
 811     printf ("Starting cost is %.0f\n", ComputeCost (Nets, T0, 5));
 812     C0 = ComputeCost (Nets, T0, T);
 813     while (1)
 814       {
 815         double Cprime;
 816         pt = createPerturbation (&Selected, T);
 817         doPerturb (&pt, false);
 818         Cprime = ComputeCost (Nets, T0, T);
 819         if (Cprime < C0)
 820           {                     /* good move! */
 821             C0 = Cprime;
 822             good_moves++;
 823             steps++;
 824           }
 825         else if ((random () / (double) RAND_MAX) <
 826                  exp (MIN (MAX (-20, (C0 - Cprime) / T), 20)))
 827           {
 828             /* not good but keep it anyway */
 829             C0 = Cprime;
 830             steps++;
 831           }
 832         else
 833           doPerturb (&pt, true);        /* undo last change */
 834         moves++;
 835         /* are we at the end of a stage? */
 836         if (good_moves >= good_move_cutoff || moves >= move_cutoff)
 837           {
 838             printf ("END OF STAGE: COST %.0f\t"
 839                     "GOOD_MOVES %d\tMOVES %d\t"
 840                     "T: %.1f\n", C0, good_moves, moves, T);
 841             /* is this the end? */
 842             if (T < 5 || good_moves < moves / CostParameter.good_ratio)
 843               break;
 844             /* nope, adjust T and continue */
 845             moves = good_moves = 0;
 846             T *= CostParameter.gamma;
 847             /* cost is T dependent, so recompute */
 848             C0 = ComputeCost (Nets, T0, T);
 849           }
 850       }
 851     changed = (steps > 0);
 852   }
 853 done:
 854   if (changed)
 855     {
 856       DeleteRats (false);
 857       AddAllRats (false, NULL);
 858       ClearAndRedrawOutput ();
 859     }
 860   FreePointerListMemory (&Selected);
 861   return (changed);
 862 }