//  The "simulated annealing" algorithm was developed by Nicholas Metropolis in 1953.
//  It was featured in "The Amateur Scientist" column by Shawn Carlson, in Scientific
//  American, March, 1997.

#include "common.h"
#include "cstring.h"
#include "cyber.h"
#include "disk.h"
#include "floating.h"
#include "anneal.h"


  inline bool boltzmann(real E, real kT) {
if (E<0.0) return 1; return (ulong)Random() < 65536L * exp(-E / kT); }

enum { alterNone,alterFlip,alterMove };

typedef struct {
short type,    //  whether itŐs a reversal, or a cut-and-paste
      start,   //  index of first node of selection
      end,     //  first node after selection
      insert;  //  first node after insertion point, if applicable
} alteration;

  static void randalteration(int NUM, alteration *A) {
if (!A) return; A->type=alterNone; if (NUM<3) return;
A->start=randzb(NUM); A->end= randzb(NUM) + 1;
if (A->start>=A->end) { if (onein(2)) A->start=0; else A->end=NUM; }
if (!range(A->start,A->end) || A->end>NUM) error("randalteration() range error #1");
A->insert=randzb(NUM + 1);
if (A->insert<A->start || A->insert>A->end) A->type=alterMove;
else A->type=alterFlip; }

  static void alterlist(annealtask *T, alteration *A) {
short i,j,a,b,c; void *TEMP;
if (A->type==alterFlip) {
  for (i=A->start,j= A->end - 1;i<j;i++,j--) {
  TEMP=T->order[i]; T->order[i]=T->order[j]; T->order[j]=TEMP; }
  return; }
if (A->type==alterMove) {
  if (A->insert<A->start) { a=A->insert; b=A->start; c=A->end; }
  else { a=A->start; b=A->end; c=A->insert; }
  i=a; j= a + c - b; while (i<b) T->space[j++]=T->order[i++];
  i=b; j= a;         while (i<c) T->space[j++]=T->order[i++];
  for (i=a;i<c;i++) T->order[i]=T->space[i]; return; }
error("alterlist() error #1"); }

  static real changeincost(annealtask *T, alteration *A) {
real Old,New; void **S,**E,**I; annealtermpp EC;
if (!T || !A || A->type==alterNone) return 0.0;
Old=0.0; New=0.0; EC=T->endcost;
S= T->order + A->start; E= T->order + A->end; I= T->order + A->insert;

if (A->start>0) Old+=T->cost(S[-1],S[0]); else if (EC) Old+=EC(S[0]);
if (A->end<T->num) Old+=T->cost(E[-1],E[0]); else if (EC) Old+=EC(E[-1]);

if (A->type==alterFlip) {
  if (A->start>0) New+=T->cost(S[-1],E[-1]); else if (EC) New+=EC(E[-1]);
  if (A->end<T->num) New+=T->cost(S[0],E[0]); else if (EC) New+=EC(S[0]); }

else if (A->type==alterMove) {
  if (A->insert>0 && A->insert<T->num) Old+=T->cost(I[-1],I[0]);
    else if (EC) { if (A->insert==0) Old+=EC(I[0]);
    else if (A->insert==T->num) Old+=EC(I[-1]); }
  if (A->start>0 && A->end<T->num) New+=T->cost(S[-1],E[0]);
    else if (EC) { if (A->start==0) New+=EC(E[0]);
    else if (A->end==T->num) New+=EC(S[-1]); }
  if (A->insert>0) New+=T->cost(I[-1],S[0]); else if (EC) New+=EC(S[0]);
  if (A->insert<T->num) New+=T->cost(E[-1],I[0]); else if (EC) New+=EC(E[-1]); }

else error("changeincost() error #1");
return New - Old; }

  inline bool annealattemp(annealtask *T) {
long i,SUCC; alteration A; real DeltaE; SUCC=0;
for (i=0;i<T->tries;i++) {
  randalteration(T->num,&A); DeltaE=changeincost(T,&A);
  if (!boltzmann(DeltaE,T->temp)) continue;
  alterlist(T,&A); SUCC++; if (SUCC>=T->oodles) return 1; }
return SUCC>0; }

  static void initanneal(annealtask *T, annealcostpp Cost, annealtermpp EndCost) {
int i; real Energy; if (!T || !Cost) return;
T->cost=Cost; T->endcost=EndCost; Energy=0.0;
for (i=0;i < T->num - 1;i++) Energy+=T->cost(T->order[i],T->order[i + 1]);
T->temp= Energy / T->num; T->tries= 100L * T->num; T->oodles= 10L * T->num; }

  void anneal(annealtask *T, annealcostpp Cost, annealtermpp EndCost) {
int i; if (!T) return; initanneal(T,Cost,EndCost);
for (i=0;i<100;i++) { if (!annealattemp(T)) return; T->temp*=0.9; } }