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