// An implementation of The Algorithm of the Gods.
//
//	Copyright 1997
//	Eric McRae	eric@elmi.com
//	Electro-Logic Machines, Inc.
//	Port Townsend, WA
//
// Original in C written by Shawn Carlson and available at
// http://web2.thesphere.com/SAS/SciAm/1997/SimAneal/SimAneal.html
//
// You may use this code as you see fit and you may redistribute it in
// whole or in part as long as the above copyright remains with any code
// you use.
// 
// The algorithm for the fairly random number generator is from my friend
// Tim Nelson.
//
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

// You may have to adjust the following for your machine/compiler
// Just make sure that the size of the variable matches the number of bits
// in the type.
typedef unsigned char u8;
typedef unsigned long u32;

// This value represents the maximum possible length of an input list.
// Once the list has been read, all other functions use the actual Item count.
#define MAXITEMS 200

// The temperature reduction factor.  The temperature is multiplied by
// this amount after each anneal pass.
#define TEMPFACTOR 0.9557F

// The list alteration type enum
enum ALTERATION
{
    REVERSE,		// reverse the endpoints of the list segment
    MOVE		// Move the segment
};

// a binary enum
enum CHOICE
{
    NO,
    YES
};

// The Item class.  This class and it's members are peculiar
// to the input type and should be modified to match your needs.
// In the case of this program, an item is just a point in an X-Y plane.
class Item
{
    public:
    				// the Item constructor
	Item( float x=0.0, float y=0.0) : xF(x), yF(y) {};

				// calculates the energy between two items
	friend float linkEnergyF( Item *, Item * );

				// display the Item values on an output stream
	friend ostream & operator << (ostream &os, Item i )
	    { return( os << i.xF << '\t' << i.yF ); }

    	const float xF, yF;	// the coordinates of the point
};


// The ItemList class handles a list of Items
class ItemList
{
    public:
	ItemList()		// ItemList Constructor
	    { fillPP = itemsAP; savedEnergyF = 10000000.0F; actual = NO; }

				// add a new Item to the list
	void addItemV( Item *iP ) { *fillPP++ = iP; }


	    			// the Anneal Function
	int annealN( void );

				// pick a random segment and insert point
	void pickSegmentV( ALTERATION );

				// alter the list as picked by pickSegmentV()
	void alterListV( ALTERATION );

				// Calc the energy change due to an alteration
	float deltaEnergyF( ALTERATION );

				// Get the energy in the list
	float getEnergyF( CHOICE =NO );

				// get the energy in the saved list
	float getSavedEnergyF( void );

				// display the list values on an output stream
	friend ostream & operator << (ostream &, ItemList & );

				// read an input stream into list
	friend istream & operator >> (istream &, ItemList & );

    private:
	CHOICE actual;		// use pseudo or actual energy
	void copyListV( void );	// for copying the best list so far
				// calculates the energy between two items
	float linkEnergyF( Item *, Item * );
	Item **startPP, **endPP, **insertPP;	// from pickSegment
	int segSizeN;				// from pickSegment
	Item *itemsAP[MAXITEMS];		// array of Item pointers
	Item *savedItemsAP[MAXITEMS];		// array of Item pointers
	float savedEnergyF;		// the lowest pesudo energy so far
	Item **fillPP;				// fill pointer for new Items
};


// statics
static int succLimitN;		// the max number of changes allowed/anneal
static int itemCountN;		// the total number of Items read in
static u32 x, y, z;		// used in the randB()
static float temperatureF;	// the annealing temperature
static ItemList list;		// The list object
static int runLimitN;		// max anneal attempt per temperature


// Functions, methods, etc.

// Function:	seedRandV
// Description:	Seeds the fairly random number generator
// Arguments:	u32 Seed, should have many bits set
// Returns:	void
void
seedRandV( u32 valL )
{
    x = valL;
    y = valL << 8;
    z = valL >> 8;
}


// The maximum value of the fairly random number generator
#define MAXRAND 255

// Function:	randB
// Description:	fairly random number generator
// Arguments:	none
// Returns:	u8 value between 0 and MAXRAND
u8
randB( void )
{
    register u32 tmpL;

    tmpL = ((x+=y+=z+=x) & 0x7fffffff) >> 19;

    return( (u8)tmpL );
}


// Function:	oracle
// Description:	decides whether to use a particular positive energy
//		modification
// Arguments:	float energy difference of modification
//		float temperature
// Returns:	CHOICE, YES = do the modification
// Notes:	The exp math function is very expensive performance wise so
//		tests are made for arguments which will always result in YES
//		or NO.  The exp function is avoided in these cases.
CHOICE
oracle( float ediF, float tempF )
{
    register float facF = ediF/(tempF*itemCountN);

    if( facF < 0.003929F )		// exp(-0.003929) = 254/255
	return( YES );

    else if( facF > 5.541264F)		// exp(-5.541264) = 1/255
	return( NO );

    else if( (int)( (exp(-facF)) * MAXRAND ) > randB() )
	return( YES );

    else
	return( NO );
}


// Function:	main
// Description:	reads an input file and writes an output file.  The input
// file contains a list of Items.  The output file gets the resulting list
// with the lowest discovered energy
int
main( int argc, char **argv )
{
    ifstream infile;
    ofstream outfile;
    int changesN;		// number of actual changes during anneal
    float energyF, bestEnergyF;

    				// check out input and output files
    switch( argc )		// switch leaves room for adding args such as 
    {				// TEMPFACTOR and others.
    case 3:
	if( strcmp( *(argv+1), *(argv+2) ) == 0 )
	{
	    cerr << *argv <<
		": input and output files cannot be the same" << endl;
	    return( -1 );
	}
	infile.open( *(argv+1) );	// try to open both input and output
	outfile.open( *(argv+2) );
	break;

    default:
	cerr << "Usage: " << *argv << " infile outfile" << endl;
	return( -2 );
    }

    if( ! infile.good() )	// if we failed to open the input file
    {
	cerr << *argv << ": can't open input file: " << *(argv+1) << endl;
	return( -3 );
    }

    if( ! outfile.good() )	// if we failed to open the output file
    {
	cerr << *argv << ": can't open output file: " << *(argv+2) << endl;
	return( -4 );
    }
    outfile.close();		// will be opened and closed later as needed


				    // set output stream format
    outfile.setf(ios::fixed,ios::floatfield);	// fixed point
    outfile.setf(ios::right,ios::adjustfield);	// right justified

				    // set stdout stream format too
    cout.setf(ios::fixed,ios::floatfield);	// fixed point
    cout.setf(ios::right,ios::adjustfield);	// right justified

    // Now read the input file and create an array of pointers to
    // each of the Items.

    infile >> list;

    		// now that we know how many items we have, initialize these
    		// globals so other functions can make fast comparisons
    succLimitN = 10 * itemCountN;
    runLimitN = 100 * itemCountN;

    infile.close();			// All done with the input file

    // OK now we've read in the Items.  Let's get on with the process.

    seedRandV(0x5a5a5a5a);		// seed the random number generator

    energyF = list.getEnergyF();	// get starting energy of list

    temperatureF = energyF/itemCountN;	// create a starting temperature

		    // Pseudo energy is displayed/used until the very end.
    //cout << "Pass\tTemp\tP-Energy\tChanges" << endl;

    bestEnergyF = 10000000.0F;		// will be replaced on the first pass

    for( register int i = 0; 1; i++ )	// go forever
    {
	changesN = list.annealN();	// anneal the list

	energyF = list.getEnergyF();	// get it's pseudo energy

				// display something to entertain the human
	//cout << i + 1 << '\t' <<
	    //temperatureF << '\t' << energyF << '\t' << changesN << endl;

	if( changesN == 0 )	// if anneal couldn't, we're done
#if 0
	    break;
#else
	{
	    temperatureF = energyF/itemCountN;
	    cout << "Pseudo Energy: " << energyF << endl;
	    energyF = list.getSavedEnergyF();	// get and display the lowest energy
	    if( energyF < bestEnergyF )
	    {
		bestEnergyF = energyF;	// if just discovered a better list
		outfile.open( *(argv+2) );
		outfile << list;
		outfile.close();
		cout << "Best Actual Energy: " << energyF << endl;
	    }
	}
	else
#endif
	temperatureF *= TEMPFACTOR;	// otherwise, reduce temperature
    }

    				// Done annealing, report on the best path
    outfile << list;		// print the best list to the output file
    outfile.close();

    energyF = list.getSavedEnergyF();	// get and display the lowest energy
    cout << "Best Actual Energy: " << energyF << endl;

    return(0);
}


// Method:	annealN
// Description:	anneals a list of Item pointers if possible
// Arguments:	none
// Returns:	number of alterations made
// Notes:
int
ItemList::annealN( void )
{
    register int i, changesN = 0;
    register ALTERATION alteration;
    register float energyDiF;

    for( i = runLimitN; i; i-- )	// attempt a bunch of changes
    {
			// randomly chose to reverse a segment or move it
	alteration = randB() & 1 ? REVERSE : MOVE;

			// randomly pick a segment to be manipulated
	list.pickSegmentV( alteration );

			// see what the energy difference would be
	energyDiF = list.deltaEnergyF( alteration );

			// should we use this alteration?
	if( energyDiF < 0 )
	{		// Always use anything that results in lower energy
	    list.alterListV( alteration );
	    changesN++;
	}
	else if( oracle( energyDiF, temperatureF ) )
	{		// go for higher energy only if oracle says so
	    list.alterListV( alteration );
	    changesN++;
	}

	if( changesN >= succLimitN )	// if it's too hot!
	    break;
    }
    return( changesN );
}


// Method:	alterListV
// Description:	asserts the previously calculated alteration
// Arguments:	ALTERATION choice to be used
// Returns:	void
// Notes:	re-arranges the Item pointer list
void
ItemList::alterListV( ALTERATION alter )
{
    Item *holdAP[MAXITEMS];	// temp storage for block of pointers
    register Item *holdP;	// temp storage for one pointer
    register Item **upPP, **downPP, **srcPP, **dstPP;	// working ptrs

    if( alter == REVERSE )
    {				// just reversing a segment in place
	upPP = startPP;
	downPP = endPP;

	while( upPP < downPP )
	{
	    holdP = *upPP;
	    *upPP++ = *downPP;
	    *downPP-- = holdP;
	}
    }

    else	// alter == MOVE, moving a segment somewhere.  For moves
    {		// up, the block is inserted above the insert point.  For
		// moves down, the block is inserted below the insert point.

		// first copy out the block we're going to move
	dstPP = holdAP;
	srcPP = startPP;
	while( srcPP <= endPP )
	    *dstPP++ = *srcPP++;

	if( insertPP > endPP )	// if insert point is above block
	{
	    dstPP = startPP;	// move from insert point down
	    srcPP = endPP + 1;
	    while( srcPP <= insertPP )
		*dstPP++ = *srcPP++;

	    srcPP = holdAP;	// then copy back the saved block
	    dstPP = insertPP - segSizeN + 1;
	    for( int i = segSizeN; i; i-- )
		*dstPP++ = *srcPP++;
	}
	else			// if the insert point is below the block
	{
	    dstPP = endPP;	// move pointers up into space left by block
	    srcPP = startPP - 1;
	    while( srcPP >= insertPP )
		*dstPP-- = *srcPP--;

	    srcPP = holdAP;	// then copy back the saved block
	    dstPP = insertPP;
	    for( int i = segSizeN; i; i-- )
		*dstPP++ = *srcPP++;
	}
    }		// end of if alteration is a MOVE
}


// Method:	deltaEnergyF
// Description:	Calculates the energy change associated with an alteration
// Arguments:	ALTERATION choice to be used
// Returns:	float energy delta
float
ItemList::deltaEnergyF( ALTERATION alter )
{
    register float costF = 0;
    register Item **lastPP = itemsAP + (itemCountN - 1);

    				// always remove energy from these links
    if( startPP > itemsAP )
	costF -= linkEnergyF( *(startPP - 1), *startPP );

    if( endPP < lastPP )
	costF -= linkEnergyF( *endPP, *(endPP + 1) );


    if( alter == MOVE )		// If we're thinking of moving a block
    {	
				// remove this energy too
	if( insertPP > itemsAP )
	    costF -= linkEnergyF( *insertPP, *(insertPP - 1) );

				// add energy associated with new situation
	costF += linkEnergyF( *insertPP, *endPP );

	if( insertPP > itemsAP )
	    costF += linkEnergyF( *(insertPP - 1), *startPP );

	if( (startPP > itemsAP ) && ( endPP < lastPP ) )
	    costF += linkEnergyF( *(startPP - 1), *(endPP + 1) );
    }

    else			// if alteration = REVERSE,
    {				// just add the endpoint changes
	if( startPP > itemsAP )
	    costF += linkEnergyF( *(startPP - 1), *endPP );

	if( endPP < itemsAP + (itemCountN  - 1) )
	    costF += linkEnergyF( *startPP, *(endPP + 1) );
    }

    return( costF );	// Note will be divided by ItemCountN in oracle
}


// Method:	pickSegmentV
// Description:	randomly picks a block of pointers and an insert point
// Arguments:	ALTERATION
// Returns:	void
// Notes:	Sets private members start, end, and insertPoint
//		The insert point is always outside the block.
//		If the alteration is a REVERSE, the insert point is not
//		updated.
void
ItemList::pickSegmentV( ALTERATION alt )
{
    register int shuffle0, shuffle1, shuffle2, hold;

    // Need to select three mutually exclusive indicies

			// first just pick a number
    shuffle0 = (randB() * (itemCountN - 1))/MAXRAND;

			// then pick another not equal to the first
    while( (shuffle1 = (randB() * (itemCountN - 1))/MAXRAND) == shuffle0 ) ;

    			// sort the first two
    if( shuffle0 > shuffle1 )
    {
	hold = shuffle1;
	shuffle1 = shuffle0;
	shuffle0 = hold;
    }

    			// If we're only reversing the segment, we do not
    			// need the insert point
    if( alt == REVERSE )
    {
	startPP = itemsAP + shuffle0;
	endPP = itemsAP + shuffle1;
    }
    else		// We're moving the segment
    {
	do		// pick a mutually exclusive third point
	    shuffle2 = (randB() * (itemCountN - 1))/MAXRAND;
	while( (shuffle0 == shuffle2) || (shuffle1 == shuffle2) );

	// Now sort these three in ascending order
	
	if( shuffle2 < shuffle0 )
	{
	    hold = shuffle0;
	    shuffle0 = shuffle2;
	    shuffle2 = shuffle1;
	    shuffle1 = hold;
	}
	else if( shuffle2 < shuffle1 )
	{
	    hold = shuffle2;
	    shuffle2 = shuffle1;
	    shuffle1 = hold;
	}

			// Randomly decide if the insert point is above
			// or below the segment.  each shuffle is random
			// so just use one of them to decide.
	if( shuffle1 > itemCountN/2 )
	{			// lets make the insert point below the block
	    insertPP = itemsAP + shuffle0;
	    startPP = itemsAP + shuffle1;
	    endPP = itemsAP + shuffle2;
	}
	else
	{			// OK, put the insert point above the block
	    startPP = itemsAP + shuffle0;
	    endPP = itemsAP + shuffle1;
	    insertPP = itemsAP + shuffle2;
	}
    }			// end of if moving a block

    segSizeN = endPP - startPP + 1;	// remember the size of this block
}


// Method:	getEnergyF
// Description:	Calculates the energy of the active or the saved Item list
// Arguments:	CHOICE  Yes if want energy from savedItemsAP
// Returns:	Total energy contained in list
// Notes:	The function declaration in the ItemList class defaults the
//		CHOICE argument to NO.
//		If the choice is NO, the function also checks to see if the
//		calculated energy is the best yet and if so, copies the active
//		list to the saved list.
//		savedEnergy is in pseudo energy units!
float
ItemList::getEnergyF( CHOICE choice )
{
    register int i;
    register float energyF = 0.0;
    register Item **listPP;

    if( choice == YES )		// if wants energy from saved list
	listPP = savedItemsAP;
    else			// Otherwise point at the working list
	listPP = itemsAP;

    for( i = 0; i < itemCountN - 1; i++ )
    {
	energyF += linkEnergyF( *listPP, *(listPP+1) );
	listPP++;
    }

    if( choice == NO )		// if this is a regular energy pass
    {				// see if this is an energy improvement
	if( energyF < savedEnergyF )	// if so
	{
	    savedEnergyF = energyF;	// remember this pseudo energy
	    copyListV();		// remember the best list
	}
    }
    return( energyF );
}


// Method:	linkEnergyF
// Description:	Calculates the energy between two Items
// Arguments:	pointer to Item
//		pointer to Item
// Returns:	energy contained in link between Items
// Notes:	if actual is YES, this function calculates the length
//		of the vector between the two items.  If actual is NO
//		someting proportional to the vector length is returned.
float
ItemList::linkEnergyF( Item * aP, Item * bP )
{
    register float dxF, dyF;

    dxF = (float)fabs(bP->xF - aP->xF);		// delta X = X2 - X1
    dyF = (float)fabs(bP->yF - aP->yF);

    if( actual )			// if need actual energy
	return( (float) sqrt( (dxF * dxF) + (dyF * dyF) ) );
    else				// else use fast approximation
	return( dxF + dyF );
}


// Method:	copyListV
// Description:	makes an internal copy of the Item list and energy.
// Arguments:	none
// Returns:	void
// Notes:	Copies the active list to the saved list
void
ItemList::copyListV( void )
{
    register Item **srcPP, **dstPP;

    srcPP = itemsAP + itemCountN - 1;
    dstPP = savedItemsAP + itemCountN - 1;

    while( srcPP >= itemsAP )
	*dstPP-- = *srcPP--;
}


// Method:	getSavedEnergyF
// Description:	reports the actual energy in the saved list
// Arguments:	none
// Returns:	float total actual energy
float
ItemList::getSavedEnergyF( void )
{
    float resultF;

    actual = YES;		// causes linkEnergy to return actual energy
    resultF = getEnergyF( YES );// get actual energy of the saved list
    actual = NO;
    return( resultF );
}


// Method:	operator <<
// Description:	prints the saved list and it's energy to an output stream
// Arguments:	ostream ref
// Returns:	void
ostream &
operator << (ostream &os, ItemList &i )
{
    Item **srcPP;

    for(srcPP = i.savedItemsAP; srcPP < i.savedItemsAP + itemCountN; srcPP++ )
	os << **srcPP << endl;

    return( os );
}


// Method:	operator >>
// Description:	Reads a list from an input stream, creating Items
//		You would seriously modify this routine if you change the
//		definition of an item.
// Arguments:	istream ref
// Returns:	void
// Notes:	Updates the global itemCountN to reflect the total number of
//		list items created.
//		Assumes we can add at least MAXITEMS Items to the list.
istream &
operator >> (istream &is, ItemList &i )
{
    for(itemCountN = 0; 1; itemCountN++ )
    {
	float x, y;

	if( itemCountN > MAXITEMS )	// > because eof not seen till bad read
	{
	    cerr << "Too many items in input" << endl;
	    return( is );		// just stop reading here
	}

	is >> x >> y;		// read two floats

	if( is.eof() )		// If that read hit EOF
	    break;			// bail out of loop

		// Otherwise, create a new Item and add it to the list
	list.addItemV( new Item(x,y) );
    }

    return( is );
}