PROGRAM pack(input,output);
{Demonstration of simulated annealing applied to the rectangle 
packing problem.  Use key 'R' to raise the temperature and 'L'
to lower.}    
         
   CONST maxn=100;
         bloat=8192;
   
   
   TYPE intarray = ARRAY [1..maxn] OF integer;
   
   
   VAR opt,flip,scale,mindim,maxdim,i,j,n,x0,y0,cost: integer;
       temp,sym,seed: integer;
       clr,xdim,ydim,xorder,yorder: intarray;
       done,improvement: boolean;
   
   
   %include 'level1.inc'

   
   FUNCTION  testsymbol: integer; EXTERN;

   
   FUNCTION randomint(n: integer): [VOLATILE] integer;
   {VOLATILE, otherwise the code produced may be distinctly non-optimal.}
   BEGIN
      seed := (seed*27913 + 1) MOD 65536;
      randomint := 1 + (seed*n) DIV 65536
   END;
   

   FUNCTION randomreal: [VOLATILE] real;
   BEGIN
      randomreal := randomint(1023)/1024
   END;
   

   PROCEDURE swap(VAR i,j: integer);
      VAR t: integer;
   BEGIN
      t:=i; i:=j; j:=t
   END;
   
   
   PROCEDURE initialise;
      VAR 
         i, j, area: integer;
         bloated_area: real;
   BEGIN
      area := 0;
      bloated_area := 0.0;
      flip := 0;
      scale := 1;
      cost := maxint;
      temp := 0;
      seed := CPUTIME MOD 65536;
      FOR i:=1 TO n DO BEGIN
         xorder[i] := i;  yorder[i] := i;
         xdim[i] := mindim + randomint(maxdim - mindim + 1) - 1;
         ydim[i] := mindim + randomint(maxdim - mindim + 1) - 1;
         area := area + xdim[i]*ydim[i];
         xdim[i] := bloat*xdim[i] + randomint(128);
         ydim[i] := bloat*ydim[i] + randomint(128);
         bloated_area := bloated_area + xdim[i]*ydim[i];
         clr[i] := (i MOD 7) + 1
      END;
      opt := 0;
      WHILE sqr(opt) < area DO opt := succ(opt);
      x0 := round(sqrt(area) + 0.5);
      y0 := x0;
      clear;
      FOR i:=n DOWNTO 2 DO BEGIN
         j:=randomint(i);
         swap(xorder[i],xorder[j]);
         j:=randomint(i);
         swap(yorder[i],yorder[j])
      END
   END;
   
   
   PROCEDURE  longest_paths(VAR x_disp, y_disp: intarray; 
                                VAR xmax, ymax: integer);
      VAR 
         inv_yorder: intarray;
         i,j,b,bb,rank,right,upper: integer;
   BEGIN
      xmax:=0; ymax:=0;
      FOR i:=1 TO n DO BEGIN
         inv_yorder[yorder[i]]:=i;
         x_disp[i]:=0; 
         y_disp[i]:=0 
      END;
      FOR i:=1 TO n DO BEGIN
         b:=xorder[i];
         rank:=inv_yorder[b];
         right := x_disp[b]+xdim[b];
         upper := y_disp[b]+ydim[b];
         FOR j:=i+1 TO n DO BEGIN
            bb:=xorder[j];
            IF rank<inv_yorder[bb] THEN BEGIN
               IF upper>y_disp[bb] THEN y_disp[bb]:=upper
            END ELSE IF right>x_disp[bb] THEN x_disp[bb]:=right
         END;
         IF right>xmax THEN xmax:=right;
         IF upper>ymax THEN ymax:=upper
      END
   END;
   
   
   PROCEDURE draw_thermometer;
   {From an original idea by M. R. King.  All rights reserved.}
      VAR i,t,x_offset,y_offset,log_temp: integer;
   BEGIN
      x_offset:=810;
      FOR i:=0 TO 1 DO BEGIN
         y_offset:=512*i;
         t:=temp;
         log_temp:=1;
         WHILE t>0 DO BEGIN t := t DIV 2; log_temp := succ(log_temp) END;
         colour(7);
         ring(x_offset, y_offset+440, 3);
         colour(0);
         fill(x_offset-3, y_offset+435, x_offset+3, y_offset+439);
         colour(1);
         disc(x_offset, y_offset+80, 8);
         colour(7);
         ring(x_offset, y_offset+80, 8);
         colour(1);
         fill(x_offset-1, y_offset+88, x_offset+2, y_offset+log_temp*11+88);
         colour(0);
         fill(x_offset-1, y_offset+log_temp*11+89, x_offset+2, y_offset+440);
         colour(7);
         line(x_offset-3, y_offset+87, x_offset-3, y_offset+439);
         line(x_offset+3, y_offset+87, x_offset+3, y_offset+439)
      END
   END;
   
   
   PROCEDURE update_screen(VAR x_disp,y_disp: intarray; xmax,ymax: integer);
      VAR i,x_offset,y_offset: integer;
   
      PROCEDURE prepare_screen(width,height: integer);
      BEGIN
         flip:=1-flip;
         WHILE (scale*width <= 200) AND (scale*height <= 200) DO
            scale:=scale*2;
         IF ((scale*width+40 > 512) OR (scale*height+40 > 512))
            AND (scale>=2) THEN scale := scale DIV 2;
         halfclear(flip);
         colour(7);
         width:=width*scale;
         height:=height*scale;
         x_offset := 448 - width DIV 2;
         y_offset := 256 - height DIV 2 + 512*flip;
         fill(x_offset-2, y_offset-2, x_offset+width+3, y_offset-1);
         fill(x_offset-2, y_offset-2, x_offset-1, y_offset+height+3);
         fill(x_offset-2, y_offset+height+2, x_offset+width+3, y_offset+height+3);
         fill(x_offset+width+2, y_offset-2, x_offset+width+3, y_offset+height+3)
      END;
   
      PROCEDURE write_box(x_ord,y_ord,w,h,c: integer);
         VAR xx,yy: integer;
      BEGIN
         xx := x_offset + scale*x_ord + 2;
         yy := y_offset + scale*y_ord + 2;
         colour(c);
         fill(xx, yy, xx + scale*w - 3, yy + scale*h - 3)
      END;
   
      PROCEDURE splat;
      BEGIN
         offset(192, flip*512)
      END;
   
   BEGIN
      prepare_screen(xmax DIV bloat,ymax DIV bloat);
      FOR i:=1 TO n DO write_box(x_disp[i] DIV bloat, y_disp[i] DIV bloat,
                                 xdim[i] DIV bloat, ydim[i] DIV bloat, clr[i]);
      splat
   END;
   
   
   PROCEDURE perturb(VAR finished: boolean);
   VAR c,newcost,r,s,t,i,j,k,xmax,ymax,change: integer;
       x_pos,y_pos: intarray;
   
      FUNCTION cost_fun: integer;
      VAR x_excess,y_excess: integer;
      BEGIN
         IF xmax>x0 THEN x_excess := xmax - x0 ELSE x_excess:=0;
         IF ymax>y0 THEN y_excess := ymax - y0 ELSE y_excess:=0;
         IF x_excess>y_excess THEN cost_fun := 2*x_excess + y_excess
         ELSE cost_fun := 2*y_excess + x_excess
      END;
   
   BEGIN
      r:=randomint(5);
      CASE r OF
         1: BEGIN s:=randomint(n); swap(xdim[s],ydim[s]) END;
         2: BEGIN s:=randomint(n-1); swap(xorder[s],xorder[s+1]) END;
         3: BEGIN s:=randomint(n-1); swap(yorder[s],yorder[s+1]) END;
         4: BEGIN
               s:=randomint(n); t:=randomint(n);
               swap(xdim[s],xdim[t]);
               swap(ydim[s],ydim[t]);
               swap(clr[s],clr[t])
            END;
         5: BEGIN
               s:=randomint(n); t:=randomint(n);
               swap(xdim[s],ydim[t]);
               swap(ydim[s],xdim[t]);
               swap(clr[s],clr[t])
            END
      END;
      longest_paths(x_pos,y_pos,xmax,ymax);
      finished := (xmax DIV bloat <= opt) AND (ymax DIV bloat <= opt);
      newcost:=cost_fun;
      IF newcost<cost THEN change:=-1
      ELSE IF newcost=cost THEN change:=0
      ELSE IF temp=0 THEN change:=2
      ELSE IF ln(randomreal)*x0*temp < (cost-newcost)*5.0 THEN change:=1
      ELSE change:=2;
      IF change<0 THEN update_screen(x_pos,y_pos,xmax,ymax);
      IF change<=1 THEN cost:=newcost
      ELSE CASE r OF
         1: swap(xdim[s],ydim[s]);
         2: swap(xorder[s],xorder[s+1]);
         3: swap(yorder[s],yorder[s+1]);
         4: BEGIN
               swap(xdim[s],xdim[t]);
               swap(ydim[s],ydim[t]);
               swap(clr[s],clr[t])
            END;
         5: BEGIN
               swap(xdim[s],ydim[t]);
               swap(ydim[s],xdim[t]);
               swap(clr[s],clr[t])
            END
      END
   END;

   
BEGIN
   writeln('                --- SIMULATED ANNEALING DEMONSTRATION ---');
   writeln;
   writeln('Depress key "R" to raise the temperature, key "L" to lower.');
   writeln('Use CTRL/Y to halt the program.');
   writeln;
   REPEAT 
      write('Number of boxes (range 2--', maxn:1, ')');
      readln(n);
      writeln
   UNTIL (n >= 2) AND (n <= maxn);
   REPEAT
      write('Minimum dimension for boxes (range 1--25)');
      readln(mindim);
      writeln
   UNTIL (mindim >= 1) AND (mindim <= 25);
   REPEAT
      write('Maximum dimension for boxes (range ', mindim:1, '--25)');
      readln(maxdim);
      writeln
   UNTIL (maxdim >= mindim) AND (maxdim <= 25);
   initialise;
   draw_thermometer;
   perturb(done);
   write('Press RETURN to start');
   readln;
   writeln;
   REPEAT
      perturb(done);
      sym:=testsymbol;
      IF (sym=ord('R')) OR (sym=ord('r')) THEN BEGIN
         IF temp=0 THEN temp:=1
         ELSE IF temp <= (maxint DIV 4) THEN temp:=temp*2;
         draw_thermometer
      END ELSE IF (sym=ord('L')) OR (sym=ord('l')) THEN BEGIN
         temp:=temp DIV 2;
         draw_thermometer
      END
   UNTIL done;
   writeln('Optimal solution obtained.');
   writeln('(No smaller square will contain the given boxes.)');
END.