(* :Title: Anneal
*)
(* :Author:
	Jeffrey A. Stern
*)
(* :Summary:

        Anneal: A Simulated Annealing Algorithm for Combinatorica
           Graphs in Mathematica

This package includes a single procedure, Anneal, which minimizes the
total length of all edges in a graph, and returns a new graph.  The
locations of all the points are the same, just exchanged to minimize
the total length of the edges connecting them.  Thus, a circular graph
will stay circular, etc.

This package depends upon having Steven Skiena's package,
Combinatorica, also loaded prior to its use.  Combinatorica2.0 is
available by anonymous ftp from sbcs.sunysb.edu in the
/pub/Combinatorica directory, and is in the Packages/DiscreteMath
subdirectory of all Mathematicas 2.0 and later. He has also written a
book on this package: "Implementing Discrete Mathematics: Combinatorics
and Graph Theory with Mathematica", through Addison- Wesley, ISBN 0-
201- 50943- 1.  (The author is greatly indebted to Skiena's work
in Combinatorica and Mathematica).

Thus, Anneal only works on graphs.  Instead of working with the
Edges of a graph directly, it first converts to an AdjacencyList (from
Combinatorica) form of the graph, in order to save time with large but
sparce graphs.  The algorithm by nature is very slow, since it is
procedural, not one of Mathematica's strengths. However, it may be
also used as a teaching tool, or as an example to modify for other
purposes, or as a model for programming compiled external code.  I
adapted the general algorithm from "Numerical Recipes in C: the art of
scientific computing", by William H. Press et. al., Cambridge U.
Press. 

Any comments, bug reports, etc. should be forwarded to:

	Jeff Stern
	Program in Social Networks
	Social Science Tower
	University of California
	Irvine, CA  92717

	jstern@aris.ss.uci.edu

	(714)856-4021
*)
(* :Context: DiscreteMath`Anneal`
*)
(* :Package Version: .111	(2/1/93 Beta Release)
*)
(* :Copyright: Copyright 1993 by Jeffrey A. Stern

This package may be copied in its entirety for nonprofit purposes only.
Sale, other than for the direct cost of the media, is prohibited.  This
copyright notice must accompany all copies.

The author, makes no representations, express or implied, with respond to this
documentation, of the software it describes and contains, including
without limitations, any implied warranties of mechantability or fitness
for a particular purpose, all of which are expressly disclaimed.  The
author shall in no event be liable for any indirect, incidental, or
consequential damages.
*)
(* :History:
	Version 0.10 by Jeffrey A. Stern, 1993.
	Version 0.11 Updated by Eric Weisstein for new Combinatorica, Oct 11, 2002.
	             Also changed iterations to be a passed parameter instead of
	             being hardwired to 100
        Version 0.111 $Film and distanceTable added by Wilberd in 2003.
*)
(* :Keywords:
	adjacency, anneal, graph, isomorphism, automorphism, network,
	embedding
*)
(* :Mathematica Version: 4.2
*)



BeginPackage["DiscreteMath`Anneal`", {"DiscreteMath`Combinatorica`"}]

Anneal::usage =
  "Anneal[g,its] exchanges vertices' locations with each other repeatedly \
  to find a global minimum for the total length of all edges in graph g.  \
  It returns a new graph with the same vertices as the original, but possibly \
  permuted.  The number of iterations may be specified as the second \
  argument.  Example: h= Anneal[g].\n
  Anneal[g,distanceTable,iterations] uses the distances in the table instead
  of those in the plane."

$Film::usage = "$Film contains the pictures generated in last call of Anneal."

Begin["`Private`"]

(* Alen computes the distance between any two points *)
(* (x1,y1), (x2, y2) where a, b, c, d are x1, x2, y1, y2 *)
(* respectively *)
Alen= Compile[{a,b,c,d}, Sqrt[(b-a)^2 + (d-c)^2] ];

Anneal[G_Graph,iterations_:100]:= Module[{npoints, nover, ad, verts,
   nlimit, p, t, i, nsucc, n, de, ans, TFACTR=0.9, tG, s,
   oldp, minp, tolerance},
 e= Edges[G];
 v= Vertices[G];
 ad=  ToAdjacencyLists[G];
 adu= ToAdjacencyLists[MakeUndirected[G]];
 verts= v;
 npoints= Length[verts];
 nover= iterations npoints +1;
 nlimit= 10 npoints -1;
 p= 0.0;
 t= 0.5;
 n1= {};
 n2= {};
 $Film={};
 tolerance= iterations - 1;

(* TOLERANCE: This is set for 99 (actually 100 because of the loop
structure below) iterations, to be a pure Simulated Annealing
algorithm.  But I often set it for 9. That is, if we should get 10
iterations in a row with the same path length, I am going to cut it
off there, and assume we've reached the global minimum, instead of
waiting for it to go all the way through 100 iterations.  This is NOT
a standard part of the pure Simulated Annealing algorithm, and
actually probably defeats the purpose of getting back out of those
local minima.  BUT, I am impatient, and computer time can be
expensive.  If you want the algorithm to work, guaranteed, through the
full 100 iterations, then just leave it set to (iterations -1)
*)
                  
(* CALCULATE INITAL TOTAL LENGTH OF ALL EDGES *)
 For[i=1, i<npoints+1, i++,
  s=Length[ad[[i]] ] + 1;
  For[j=1, j<s, j++,
   p+= Alen[verts[[i,1]], verts[[ad[[i,j]],1]],
            verts[[i,2]], verts[[ad[[i,j]],2]] ];
  ]; (* FOR j *)
 ]; (* FOR i *)
 p/= 2;
 Print["Initial Path Length"];
 Print[ToString[p]];
 oldp=p; minp=p;

Print["Iter   Length    #accepted"];
Print["====   =======   ========="];

(* TRY UP TO 100 TEMPERATURE STEPS *)
 For[j=1, j<101, j++,
  nsucc=0;
  For[k=1, k<nover, k++,
   n1=0; n2=0;

(* CHOOSE 2 DIFFERENT PTS TO EXCHANGE *)
   While[(n1==n2),
    n1= Random[Integer, {1, npoints}];
    n2= Random[Integer, {1, npoints}];
   ]; (* WHILE *)

(************ START OF REVCST MODULE ************************)
(* CALCULATE THE COST OF REVERSING THESE TWO POINTS *)
   AX= adu[[n1]];  (* n1'S NEIGHBORHOOD *)
   AY= adu[[n2]];  (* n2'S NEIGHBORHOOD *)
   LAX= Length[AX]+1; (* # of points n1 is connected to, +1 *)
   LAY= Length[AY]+1; (* # of points n2 is connected to, +1 *)
   de= 0.0;           (* de is the change in energy *)
(*   Print["AX  " <> ToString[AX]];
   Print["AY  " <> ToString[AY]];
   Print["LAX " <> ToString[LAX]];
   Print["LAY " <> ToString[LAY]]; *)

(* CALC GAIN (REDUCTION IN LENGTH) FROM DISCONNECTING n1 *)
(* FROM AX *)
   For[i=1, i<LAX, i++,
    de-= Alen[verts[[n1,1]], verts[[AX[[i]],1]],
              verts[[n1,2]], verts[[AX[[i]],2]] ];
   ];
(*   Print["de:"]; *)
(*   Print[ToString[de]]; *)
(* NOW DO THE SAME FOR n2 *)
   For[i=1, i<LAY, i++,
    If[(AY[[i]] != n1),
     de-= Alen[verts[[n2,1]], verts[[AY[[i]],1]],
               verts[[n2,2]], verts[[AY[[i]],2]] ];
    ];
   ];
   (* Print[ToString[de]]; *)
(* CALCULATE LOSS (INCREASE IN LENGTH) FROM CONNECTING n1 *)
(* UP TO n2'S NEIGHBORHOOD *)
   For[i=1, i<LAY, i++,
    If[(AY[[i]]!= n1),
     de += Alen[verts[[n1,1]], verts[[AY[[i]],1]],
                verts[[n1,2]], verts[[AY[[i]],2]] ],
    (* else *)
     de += Alen[verts[[n1,1]], verts[[n2,1]],
                verts[[n1,2]], verts[[n2,2]] ];
    ];
   ];
   (* Print[ToString[de]]; *)
(* NOW CALCULATE LOSS (INCREASE IN LENGTH) FROM CONNECTING n2 *)
(* UP TO n1's NEIGHBORHOOD *)
   For[i=1, i<LAX, i++,
    If[(AX[[i]] != n2),
     de+= Alen[verts[[n2,1]], verts[[AX[[i]],1]],
               verts[[n2,2]], verts[[AX[[i]],2]] ];
    ];
   ];
(*   Print[ToString[de]]; *)
(************ END OF REVCST MODULE ****************************)
(************ START OF METROP MODULE **************************)
   If[((de<0) || (Random[] < Exp[(-de/t)])),
    ans=True, ans=False;
   ];
(************ END OF METROP MODULE ****************************)
   If[ans,
    ++nsucc;
    p+= de;
(************ START OF REVERSEPOINTS MODULE *******************)
    hv= verts[[n1]];
    verts[[n1]]= verts[[n2]];
    verts[[n2]]= hv;
(************ END OF REVERSEPOINTS MODULE *********************)
   ]; (* If *)
   If[nsucc > nlimit, Break[] ];
  ]; (* FOR k *)
  Print["   " <>
        ToString[j] <> "  " <>
        ToString[PaddedForm[p,{6,4}]] <> "     " <>
        ToString[nsucc]];
	$Film=Append[$Film,
	ShowGraph[Graph[List/@e,List/@verts],
	     PlotLabel->ToString[PaddedForm[p,{6,4}]],
	     DisplayFunction->Identity]];
	ShowGraph[Graph[List/@e,List/@verts]];

  t *= TFACTR;
  If[p==minp, ++Shorter,
   If[p<minp, Shorter=1; minp=p, Shorter=0]];
  If[((nsucc==0) || (Shorter > tolerance)), Break[] ];
 ]; (* FOR j *)
 
(* MAKE NEW GRAPH WITH NEW VERTICES AND SEND BACK *)
 tG= Graph[List/@e,List/@verts];
 Print["Final Path Length"];
 Print[ToString[p]];
 Print["DONE."];
 $Film = $Film /. Identity -> $DisplayFunction;
 tG]; (* ANNEAL *)

Anneal[g_Graph,distanceTable_List,iterations_]:=
Block[{Alen,auxGraph},
  nverts=Length[Vertices[g]];
  auxGraph=ChangeVertices[g,Table[{i,nverts Random[]}, {i, nverts}]];
  Alen[a_, b_, c_, d_] := distanceTable[[a, b]];
  auxGraph=Anneal[auxGraph,iterations];
  auxGraph=ChangeVertices[auxGraph,
            Table[Vertices[g][[Vertices[auxGraph][[i,1]]]],{i, nverts}]]
       ];

End[]

EndPackage[]

      
