/*	Steven Skiena
	Vehicle routing problem
	June 1, 2000

	This program assigns maintainence jobs among operators so as to maximize
	overall efficiency.  To compile, use:

		gcc -lm -O2 routing.c

	The basic algorithm is:

	(1) Distribute all jobs to operators based on the distance to the
	operator's focal pointer and current load. 

	(2) Use local optimization / simulated annealing to
		swap and shift jobs so as to create additional slack.

	Each operator's schedule is represented by an ordered sequence of jobs,
	which defines a computed schedule of where they are at each time.
	The computed schedule currently incorpates start position, 
	job processing and travel time.  It does not yet incorporate
	the desired end position, slack periods, and lunch.
 
	The current evaluation function for the "badness" of the schedule weighs
	both travel costs and "trouble", the latter being function of both the
	number of jobs missing desired windows, as well as by how much they miss
	the window.

	The simulated annealing moves supported allow us to:
		(1) Move a job j1 from operator A to operator B.
		(2) Swap jobs j1 and j2 for single operator A.
		(3) Swap job j1 from operator A with job j2 from operator B.

	Further improving the performance of the program can be done in several ways:
		(1) Tune the constants which regulate the number of iterations and
			the distribution of moves selected to try. 
		(2) Incorporate more powerful moves, which identify the most
			problematic jobs and the most underutilized operators.
		(3) Exhaustively optimize the final schedule for each operator
			by trying all permutations
	

	Simplifications:

	This implementation simplifies the general problem in the following
	ways which enabled me to produce a working prototype within the prescribed
	time period.

	Simplify Input Files -- to make the files I need easier to parse, I have
		deleted the first line from each file, and eliminated fields which
		I do not need.  I have also eliminated the blank from all postcodes,
		and "" surrounding fields.  These files were partially created by
		hand and partially by the included awk scripts.

	End positions -- Since not all operators have endpositions, I 
		ignored this constraint.  

	Travel times -- This is taken as only being a function of the distance,
		ignoring the river constraints in postcodesplit.csv

	Dates -- this program ignores dates, worrying only about time windows.

	Productivity and Skill sets -- this implementation ignores productivity
		and skill set differences between operators.  This can be fairly
		easily incorporated later.

	Teams -- This implementation ignores the need for each team member
		to be at the job at the same time.  We will ensure that enough
		people work at the job as needed within the time window, but overlapping
		intervals will not be assured.  This can be handled better by
		post-processing the schedules associated with each team, or making the
		time window for the shared jobs to be so tight that all are there
		together.

*/

#include <stdio.h>

#define MAXOPERATIVES		155	/* Maximum number of operators */
#define MAXJOBS			2000	/* Maximum number of jobs in the system */
#define MAXPOSTCODES		1000	/* Maximum number of postal codes */
#define MAXWORKTYPES		100	/* Maximum number of distinct work types */
#define MAXTIMEBANDS		10	/* Maximum number of time bands */

#define METERS_PER_MINUTE	750	/* Generic operator travel speed */
#define MAX_JOBS_PER_OPERATIVE	9	/* How many jobs to one guy? */
#define START_OF_DAY		(8*60)	/* minute count of when each day starts */
#define NCLOSE			50	/* Size of list of nearby operatives */
#define TRAVEL_WEIGHT		5	/* Factor to combine travel/trouble penalty */

#define IDLENGTH		30	/* Length of short strings like postcodes */
#define NULL_ID			-1	/* blank code for non-existant ID */
#define LARGE_INT		1000000	/* large enough integer to be infinity */

typedef struct timeband_struct timeband_struct; /* eliminate extra struct */
struct timeband_struct {
	int id;				/* id number of timeband */
	int start;			/* start time of band */
	int length;			/* length of band */
	char descr[IDLENGTH];		/* description string */
};

timeband_struct timeband[MAXTIMEBANDS+1];	/* array of timebands */
int ntimebands;					/* number of active time bands */

typedef struct	postcode_struct postcode_struct;	/* eliminate extra struct */
struct postcode_struct {
	char id[IDLENGTH];		/* code string */
	int x, y;			/* x,y coordinates for the location */
};

postcode_struct postcode[MAXPOSTCODES];	/* array of postcodes */
int npostcodes;				/* number of active postcodes */

typedef struct schedule_struct schedule_struct; /* eliminate extra struct */
struct schedule_struct {
	int id;				/* id number of job */
	int starttime;			/* time of starting job, after travel */
	int endtime;			/* completion time of job */
	int traveltime;			/* minutes of travel to this job */
	int trouble;			/* how early/late or badly scheduled is job? */
};

typedef struct operative_struct operative_struct; /* eliminate extra struct */
struct operative_struct {
	int id;				/* id number for operative -- note not index */
	char name[IDLENGTH];		/* operator name */
	int startx, starty;		/* x,y coordinates for the start position */
	int focusx, focusy;		/* x,y coordinates for the focus location */
	int lunchlength;		/* length of lunch period in minutes */
	int lunchperiod;		/* point in day when lunch should begin */
	int jobcount;  			/* number of jobs assigned to operative */
	int jobs[MAX_JOBS_PER_OPERATIVE];    /* array of job indicies */
	schedule_struct schedule[MAX_JOBS_PER_OPERATIVE]; /* schedule of these jobs*/
};

operative_struct operative[MAXOPERATIVES];	/* array of operatives */
int noperatives;				/* number of active operatives */

typedef struct worklist_struct worklist_struct; /* eliminate extra struct */
struct worklist_struct {
	int id;				/* id number of job -- note not index */
	char postcode[IDLENGTH];	/* postcode location -- used to find x,y */
	int x, y;			/* point location of worksite */
	int jobtype;			/* id of jobtype */
	int window;			/* time window to do the job */
	int service;			/* service level -- not currently used */
	int aux;			/* auxilliary type */
	int local[NCLOSE];		/* the closest workers */
};

worklist_struct worklist[MAXJOBS];	/* array of jobs */
int nworklist;				/* number of active jobs */

typedef struct worktype_struct worktype_struct; /* eliminate extra struct */

struct worktype_struct {
	int id;				/* id number of job -- seems to be the index */
	char descr[IDLENGTH];		/* text string describing job type */
	int length;			/* number of minutes til completion */
	int skilled, unskilled;		/* number of skilled/unskilled workers needed */
};

worktype_struct worktype[MAXWORKTYPES];	/* array of work types */
int nworktype;				/* number of active worktypes */

typedef enum { FALSE=0, TRUE=1 } bool;	/* boolean data type */



/*************************************************************************/
/**************** Simulated Annealing Constants **************************/
/*************************************************************************/

#define REPEAT_COUNT		1	/* how many solution attempts do you
					   want?  More than 1 enables you to
					   eyeball the output and pick the
					   best.  If you are getting stuck in
					   local optima, this good to try. */

#define	INITIAL_TEMPERATURE	1	/* start temperature -- probably
					   leave intact */

#define COOLING_STEPS		500	/* how many times do we cool -- make
					   higher to improve quality, lower to
					   speed the program up.  Move in
					   tandem with the COOLING_FRACTION */

#define COOLING_FRACTION	0.97	/* how much to cool each time -- make
					   higher to improve quality, lower to
					   speed the program up. */

#define STEPS_PER_TEMP		500	/* lower makes it faster, higher makes
					   it potentially better. */

#define E			2.718	/* number e -- probably leave intact*/

#define K			20.5	/* problem specific Boltzman's constant
					   May have to adjust if your global
					   value function changes the sizes
					   of the numbers it produces.  It is
					   important that jumps seem random at
					   the begining of the run, and rare
					   at the end of a run, and this is
					   a knob to tweak that. */

#define RAND_MAX	32857		/* large integer for rand number generation */

#define TRACE_OUTPUT	TRUE      	/* print the swaps as they happen */


/*	Make these probabilities sum to 1! 	*/
#define PROB_SINGLE_SWAP	0.25	/* probability of selecting a single operator swap */
#define PROB_DOUBLE_SWAP	0.25	/* probability of selecting a two operator swap */
#define PROB_SIMPLE_MOVE	0.25	/* probability of selecting a simple job move */
#define PROB_MACHO_MOVE		0.25	/* probability of selecting a macho job move */

/**************************************************************************/

main()
{
	int i;				/* counter */
	read_timebands();
	read_possibledates();
	read_postcodes();
	read_operatives();
	read_worklist();
	read_worktype();

	find_closest_operatives();
	assign_operatives_to_jobs();
	/* print_operative_assignments(); */

	schedule_operatives();

	for (i=1; i<=REPEAT_COUNT; i++) {
		anneal();
		print_operative_schedules();
	}
}

/**************************************************************************/
/*	Parsing Input Files						  */
/**************************************************************************/

/*	The files we need to read to do the job are timebands.csv operatives.csv
	worklist.csv postcodes.csv worktypes.csv -- each preprocessed into
	the new format described in the comments below
*/

/* The timebands.sim file format has one line per timeband, each with
	four fields: (1) id, (2) start time of band, (3) length of band,
	(4) description.
*/

read_timebands()
{
	FILE *fp;                       /* file to open */
	int i;                          /* counter */

	fp = fopen("timebands.sim","r");

	i=1;
	while (fscanf(fp,"%d %d %d %s\n", &timeband[i].id, &timeband[i].start,
				&timeband[i].length,&timeband[i].descr) != EOF) {
		if (timeband[i].id != i)
			printf("Error: Timeband file out of order\n");
		i++;
	}

	ntimebands = i-1;

	/*print_timebands();*/
}

print_timebands()
{
	int i;                           /* counter */

	printf("timebands\n");
	for (i=1; i<=ntimebands; i++) {
		printf("%d %d %d %s\n",timeband[i].id, timeband[i].start,
				timeband[i].length, timeband[i].descr);
   	}

   	printf("\n");
}

/*	The possibledates.sim file format has one line per job, each with two
	fields: (1) the job id, (2) the time band.

	The timeband will be associated with the given entry in worklist.
*/

read_possibledates()
{
	FILE *fp;                       /* file to open */
	int id, band;			/* placeholders */

	fp = fopen("possibledates.sim","r");

	while (fscanf(fp,"%d %d\n", &id, &band) != EOF) {
		worklist[id].window = band;
	}
	fclose(fp);
}


/*	The postcodes.sim file format has one line per postcode, each with three
	fields: (1) the postcode, (2) the x-coordinate, and (3) the y-coordinate.
*/

read_postcodes()
{
	FILE *fp;                       /* file to open */
	int i;				/* counter */

	fp = fopen("postcodes.sim","r");

	i=0;
	while (fscanf(fp,"%s %d %d\n", &postcode[i].id,
			&postcode[i].x, &postcode[i].y) != EOF) {
		i++;
	}		

	npostcodes = i;

	/*print_postcodes();*/

	fclose(fp);
}

print_postcodes()
{
	int i;				/* counter */

	printf("postcodes\n");
	for (i=0; i<npostcodes; i++) {
		printf("%s	%d	%d\n",postcode[i].id,
				postcode[i].x, postcode[i].y);
	}

	printf("\n");
}

/*	Find the index of the post code corresponding to the give string */

find_postcode(s)
char *s;				/* postcode in question */
{
	int i;				/* counter */

	for (i=0; i<npostcodes; i++) {
		if (strcmp(s,postcode[i].id) == 0) return(i);
	}

	printf("Warning: postcode %s not found\n", s);
	return(NULL_ID);
}

/*	The worklist.sim file format has one line per job, each with five fields;
	(1) the job id, (2) the location postal code, (3) work-type-id, 
	(4) the service level (presumed to be the time period to serve,
	(5) the associated work ID (currently ignored)

	The worklist will be stored such that the index matches the id.
	This means that undefined jobs will be represented by blank
	records with null ids.
*/

read_worklist()
{
	FILE *fp;			/* file to open */
	int i;				/* counter */
	int id;				/* id of current line */
	int postindex;			/* index of postcode */

	for (i=0; i<MAXJOBS; i++) worklist[i].id = NULL_ID;

	fp = fopen("worklist.sim","r");

	i=0;
	while (fscanf(fp,"%d",&id) != EOF) {
		if (worklist[id].id != EOF)
		   printf("Error: worklist %d assigned more than once.\n",id);

		worklist[id].id = id;

		fscanf(fp,"%s %d %d %d\n", &worklist[id].postcode,
			&worklist[id].jobtype,&worklist[id].service,&worklist[id].aux);
		postindex = find_postcode(worklist[id].postcode);
		worklist[id].x = postcode[postindex].x;
		worklist[id].y = postcode[postindex].y;

		i++;
	}		

	nworklist = i;

	/*print_worklist();*/

	fclose(fp);
}

print_worklist()
{
	int i;				/* counter */

	printf("worklist\n");
	for (i=0; i<MAXJOBS; i++) {
		if (worklist[i].id != NULL_ID)
			printf("%d %s %d %d %d %d %d %d\n", worklist[i].id,
				worklist[i].postcode, worklist[i].x,
				worklist[i].y, worklist[i].jobtype,
				worklist[i].window, worklist[i].service,
				worklist[i].aux);
	}

	printf("\n");
}

/*	The operatives.sim file format has one line per worker, each with 8 fields;
	(1) operater id, (2) operator name, (3-4) start x- and y-coordinates,
	(5-6) focus x-y coorindates, (7) lunch length, (8) lunch after time-band

	The operatives will be stored such that the index matches the id.
	This means that undefined operatives will be represented by blank
	records with null ids.
*/

read_operatives()
{
	FILE *fp;			/* file to open */
	int i;				/* counter */
	int id;				/* just read id */

	for (i=0; i<MAXOPERATIVES; i++) operative[i].id = NULL_ID;

	fp = fopen("operatives.sim","r");

	i=0;
	while (fscanf(fp,"%d",&id) != EOF) {
		if (operative[id].id != EOF)
			printf("Error: operative %d assigned more than once.\n",id);

		operative[id].id = id;
		fscanf(fp,"%s %d %d %d %d %d %d\n",
			&operative[id].name, &operative[id].startx,
			&operative[id].starty, &operative[id].focusx,
			&operative[id].focusy,&operative[id].lunchlength,
			&operative[id].lunchperiod);
		i++;
		operative[i].jobcount = 0;
	}

	noperatives = i;
	/*print_operatives();*/

	fclose(fp);
}

print_operatives()
{
	int i;			/* counter */

	printf("operatives\n");
	for (i=0; i<MAXWORKTYPES; i++) 
		if (operative[i].id != NULL_ID)
			printf("%d %s %d %d %d %d %d %d\n", operative[i].id,
				operative[i].name, operative[i].startx,
				operative[i].starty, operative[i].focusx,
				operative[i].focusy, operative[i].lunchlength,
				operative[i].lunchperiod);

	printf("\n");
}


/*	The worktypes.sim file format has one line per worktype, each with 5 fields;
	(1) id, (2) description, (3) job length, (4) skilled workers needed,
	(5) unskilled workers needed.

	The worktypes will be stored such that the index matches the id.
	This means that undefined worktypes will be represented by blank
	records with null ids.
*/

read_worktype()
{
	FILE *fp; 			/* file to open */
	int i;				/* counter */
	int id;				/* id code */

	for (i=0; i<MAXWORKTYPES; i++) worktype[i].id = NULL_ID;

	fp = fopen("worktypes.sim","r");

	i=0;
	while (fscanf(fp,"%d",&id) != EOF) {
		if (worktype[id].id != NULL_ID) 
			printf("Error: worktype %d assigned more than once.\n",id);

		worktype[id].id = id;
		fscanf(fp,"%s %d %d %d\n", &worktype[id].descr,
			&worktype[id].length, &worktype[id].skilled,
			&worktype[id].unskilled);
		i++;
	}		

	nworktype = i;

	/*print_worktypes();*/

	fclose(fp);
}

print_worktypes()
{
	int i;				/* counter */

	printf("worktypes\n");
	for (i=0; i<MAXWORKTYPES; i++) {
		if (worktype[i].id != NULL_ID)
			printf("%d %s %d %d %d\n", worktype[i].id, worktype[i].descr,
		   worktype[i].length, worktype[i].skilled,
			worktype[i].unskilled);
	}

	printf("\n");
}


/***********************************************************************/
/*		UTILITY FUNCTIONS 				       */
/***********************************************************************/

/*	Compute the time between two points.  The distances are apparently
	in meters, with a constant travel speed of METERS_PER_MINUTE
*/

travel_time(x1,y1,x2,y2)
int x1,y1;				/* coordinates of starting point */
int x2,y2;				/* coordinates of ending point */
{
	double dx, dy;			/* delta x and delta y */
	double dist;			/* distance */
	double sqrt();

	dx = (double) (x1-x2);
	dy = (double) (y1-y2);
	dist = sqrt(dx*dx + dy*dy);
	/*printf("dist = %f \n",dist);*/
	return((int) dist/METERS_PER_MINUTE );
}

/*	Travel time between operative and job */

op_job_dist(op,job)
int op;					/* index of operator */
int job;				/* index of job */
{
	return( travel_time( operative[op].focusx, operative[op].focusy,
			worklist[job].x, worklist[job].y ) );
}


/* Select and identify the possible local operatives for each job.
	Use a simple, insertion sort-like strategy: fill the local operator
	array with the first operatives, then repeatedly insert new operatives
	to order the closest by distance from the focus.
*/

find_closest_operatives()
{
	int job, op;			/* job and operator counters */
	int i;				/* counter */
	int cnt;			/* enough operatives in close list? */

	for (job=0; job<MAXJOBS; job++) {
	   if (worklist[job].id != NULL_ID) {
		op = 0;
		cnt = 0;
		while (op < MAXOPERATIVES) {
		   if (operative[op].id != NULL_ID) {
			if (cnt < NCLOSE-1) {
				insert_local_operative(job,op,cnt);
				cnt ++;
			}
			else {
				if (op_job_dist(op,job) <
					op_job_dist(worklist[job].local[NCLOSE-1],job) )
					insert_local_operative(job,op,NCLOSE-1);
			}
		   }
		   op ++;
		}

		/* print the list of closest operatives */
/*
		printf("Closest operatives for job %d at (%d,%d): \n",job,
				worklist[job].x, worklist[job].y);
		for (i=0; i<NCLOSE; i++)
		  printf("%d. operative %d at distance %d\n",i,worklist[job].local[i],
			op_job_dist(worklist[job].local[i],job));
		printf("\n\n");
*/
	   }
	}
}

/*	Insert the new operative in the list and let him slither to his correct
	position in sorted order.
*/

insert_local_operative(job,op,index)
int job;				/* job requiring operatives */
int op;					/* candidate operative */
int index;				/* position to put him */
{
	int i;				/* counter */

	worklist[job].local[index] = op;
	i = index;
	while ((i>0) && (op_job_dist(worklist[job].local[i],job) <
		op_job_dist(worklist[job].local[i-1],job)) ) {
		swap_ints(&worklist[job].local[i], &worklist[job].local[i-1]);
		i--;
	}
}

/*	Swap two integers pointed at by the arguments */

swap_ints(x,y)
int *x, *y;
{
	int tmp;			/* placeholder */

	tmp = *x;
	*x = *y;
	*y = tmp;
}

/**********************************************************************/

/*  Here we assign each job to an initial operative.  We cruise down the
    sorted list of closest operatives to each job, and assign it to the
    closest non-overloaded people.
*/

assign_operatives_to_jobs()
{
	int job;			/* counter for job to assign */
	int i, j;			/* counters */
	int jobtype;			/* type of work to be performed */
	int *freeops;			/* pointer to list of potential operators */ 
	int nworkers;			/* how many total workers, skilled and unskilled? */
	int trynext;			/* index of worker to try next */

	for (job=0; job<MAXJOBS; job++) {
	   if (worklist[job].id != NULL_ID) {
		jobtype = worklist[job].jobtype;
		/*printf("id=%d jobtype = %d\n",worklist[job].id,jobtype);*/
		freeops = &worklist[job].local[0];
		nworkers = worktype[jobtype].skilled + worktype[jobtype].unskilled;

		if (nworkers > 1)
			printf("Multiple workers (%d) for job %d\n", nworkers, job);

		trynext = 0;

		/* assign the skilled and unskilled workers */
		for (i=1; i<=nworkers; i++) {
			while ((trynext < NCLOSE) && loaded_enough(freeops[trynext]))
				trynext ++;

			if (trynext >= NCLOSE) 
				printf("Error: out of local workers for job %d\n",job);
			else {
				assign_operative_to_job(freeops[trynext],job);
				trynext ++;
			}
		}
	   }
	}
}

/*	Assign a particular operative to a particular job */

assign_operative_to_job(op,job)
int op;					/* operative */
int job;				/* job assigned to them */
{
	if (operative[op].jobcount > MAX_JOBS_PER_OPERATIVE)
		printf("Error: operative %d has too many jobs before %d\n",op,job);
	else {
		operative[op].jobs[ operative[op].jobcount ] = job;
		operative[op].jobcount ++;
	}
}

/*	Test if a particular worker is overloaded with jobs */

loaded_enough(op)
int op;					/* operative in question */
{
	return (operative[op].jobcount >= MAX_JOBS_PER_OPERATIVE-1);
}

overloaded(op)
int op;					/* operative in question */
{
	return (operative[op].jobcount >= MAX_JOBS_PER_OPERATIVE);
}


print_operative_assignments()
{
	int op,j;			/* counters */

	for (op=0; op<MAXOPERATIVES; op++) {
	   if (operative[op].id != NULL_ID) {
		printf("Operative %i -- id=%d:\n",op,operative[op].id);
		for (j=0; j<operative[op].jobcount; j++)
			printf("%d ",operative[op].jobs[j]);
		printf("\n\n");
	   }
	}
}

/*************************************************************************/

/*	These routines maintain the operatives schedules, always working from
	the ordered list of jobs for each operative.  From this ordering we
	fill in the start, end, and travel times.
*/

schedule_operatives()
{
	int op;				/* operative counter */

	for (op=0; op<MAXOPERATIVES; op++) 
		if (operative[op].id != NULL_ID)
			fill_in_schedule_for_operative(op);
}

/*	Fill out a job list into a schedule, preserving the order of jobs
	on the list.  We fill in from the start of the day.
	We incorporate no delays to wait for windows to start and ignore lunch.
*/

fill_in_schedule_for_operative(op)
int op;					/* operative whose schedule is being filled */
{
	int i;				/* counter */
	int jobid;			/* placeholder for current job */
	int previous_jobid;		/* placeholder for previous job */
	int jobtype;			/* the type of work performed */
	schedule_struct *pt;		/* pointer to schedule */

	pt = &(operative[op].schedule[0]);

	for (i=0; i<MAX_JOBS_PER_OPERATIVE; i++)
			pt[i].id = NULL_ID;

	/* fill in initial job */
	if (operative[op].jobcount >= 1) {
		jobid = operative[op].jobs[0];
		jobtype = worklist[jobid].jobtype;
		pt[0].id = jobid;
		pt[0].traveltime = travel_time(operative[op].startx,operative[op].starty,
				worklist[jobid].x, worklist[jobid].y);
		pt[0].starttime = START_OF_DAY + pt[0].traveltime;
		pt[0].endtime = pt[0].starttime + worktype[jobtype].length;
		previous_jobid = jobid;
	}
		
	/* fill in subsequent jobs */
	for (i=1; i<operative[op].jobcount; i++) {
		jobid = operative[op].jobs[i];
		jobtype = worklist[jobid].jobtype;
		pt[i].id = jobid;
		pt[i].traveltime = travel_time(worklist[jobid].x, worklist[jobid].y,
	            worklist[previous_jobid].x, worklist[previous_jobid].y);
		pt[i].starttime = pt[i-1].endtime + pt[i].traveltime;
		pt[i].endtime = pt[i].starttime + worktype[jobtype].length;
		pt[i].trouble = in_band(pt[i].starttime, worklist[jobid].window);
		previous_jobid = jobid;
	}
}

/*	Is the job in the right timeband?  Return 0 if in band, a positive
	number if it is too late, and a negative number if it is too early.
*/

int in_band(start,band)
int start;				/* start time of job */
int band;				/* time band id */
{
	int win_start, win_end;		/* start and end times for window */

	win_start = timeband[band].start;
	win_end = timeband[band].start + timeband[band].length;

	if (start < win_start) return(win_start - start);
	else if (start > win_end) return(start - win_end);
	else return(0);
}

/*	Add up the travel time incurred on a particular schedule. */

schedule_travel(op)
int op;					/* operative */
{
	int tot;			/* total travel */
	int i;				/* counter */

	tot = 0;
	for (i=0; i<operative[op].jobcount; i++)
		tot += operative[op].schedule[i].traveltime;

	return(tot);
}

/* Add up the travel trouble incurred on a particular schedule. */

schedule_trouble(op)
int op;					/* operative */
{
	int tot;			/* total trouble */
	int i;				/* counter */

	tot = 0;
	for (i=0; i<operative[op].jobcount; i++)
	      tot += abs( operative[op].schedule[i].trouble );

	return(tot);
}

/* Add up the total badness of a particular schedule */

schedule_badness(op)
int op;					/* operative */
{
	return( schedule_trouble(op) +
		(TRAVEL_WEIGHT * schedule_travel(op)) );
}

/* Add the total, global badness over all operatives */

global_value()
{
	int tot;			/* total badness */
	int op;				/* operative in question */

	tot = 0;

	for (op=0; op<MAXOPERATIVES; op++)
		if (operative[op].id != NULL_ID)
			tot += schedule_badness(op);

	return(tot);
}

/* Add the total, travel badness over all operatives */

global_travel()
{
        int tot;			/* total badness */
        int op;				/* operative in question */

        tot = 0;

        for (op=0; op<MAXOPERATIVES; op++)
                        if (operative[op].id != NULL_ID)
                                tot += schedule_travel(op);

        return(tot);
}

global_trouble()
{
        int tot;			/* total badness */
        int op;				/* operative in question */

        tot = 0;

        for (op=0; op<MAXOPERATIVES; op++)
                        if (operative[op].id != NULL_ID)
                                tot += schedule_trouble(op);

        return(tot);
}


print_schedule_for_operative(op)
int op;					/* operative in question */
{
	int i;				/* counter */
	schedule_struct *pt;		/* shorthand for schedule record */

	printf("\nOperative %d -- badness=%d travel=%d trouble=%d\n",op,
		schedule_badness(op), schedule_travel(op), schedule_trouble(op));
	printf("id	  start  end	travel trouble\n");

	pt = &(operative[op].schedule[0]);

	for (i=0; i<MAX_JOBS_PER_OPERATIVE; i++) {
		if (pt[i].id != NULL_ID)
			printf("%d	%d	%d	%d		%d\n",
				pt[i].id, pt[i].starttime, pt[i].endtime,
				pt[i].traveltime, pt[i].trouble);
	}

	printf("\n");
}

print_operative_schedules()
{
   int op;                    /* operative counter */

	   for (op=0; op<MAXOPERATIVES; op++) 
		      if (operative[op].id != NULL_ID)
					print_schedule_for_operative(op);
}


/***********************************************************************/

/*	These routines implement simulated annealing.  Pairs of elements
	of the same type will be swapped at random, and the new arrangment
	accepted either if (1) it is an improvement, or (2) the penalty is
	less than a random number, which is a function of the temperature
	of the system.

	This is a fairly general simulated annealing routine.  The task for
	the user is to correctly implement the swap and function evaluations.

	We are seeking to minimize the current_value.
*/

anneal()
{
	double temperature;		/* the current system temp */
	int i,j;			/* counters */
	int start_value;		/* initial value at start of loop */
	int current_value;		/* value of current state */
	int current_local_value,
	    new_local_value;		/* local value before/after swap */
	int delta;			/* change in value from current to new */
	double merit, flip;		/* hold swap accept conditions*/
	double exponent;		/* exponent for energy funct*/
	bool swapflag;			/* do I swap jobs or move them? */
	bool machoflag;			/* are we doing a macho move? */
	int rand_op1,rand_op2;		/* operators to permute */
	int rand_job1, rand_job2;	/* pair of jobs to permute */
	double r;			/* hold random integer */
	double random_float(), pow();

	temperature = INITIAL_TEMPERATURE;

	current_value = global_value();

	for (i=1; i<=COOLING_STEPS; i++) {
		temperature *= COOLING_FRACTION;

		start_value = current_value;		
		printf("current_value = %d\n", current_value);
		if (TRACE_OUTPUT) {
			printf("global_travel= %d  global_trouble=%d\n",
			global_travel(), global_trouble());
		}
		
		for (j=1; j<=STEPS_PER_TEMP; j++) {

			swapflag = TRUE;
			machoflag = FALSE;
			r = random_float(0,1);
			if (r <= PROB_SINGLE_SWAP) 
			   single_op_swap(&rand_op1, &rand_job1, &rand_op2, &rand_job2);
			else {
			   r -= PROB_SINGLE_SWAP;
			   if (r <= PROB_DOUBLE_SWAP)
				double_op_swap(&rand_op1, &rand_job1, &rand_op2, &rand_job2);
			   else {
				swapflag = FALSE;		/* move a job instead of swap it */
				r -= PROB_DOUBLE_SWAP;
				if (r <= PROB_SIMPLE_MOVE)
					double_op_move(&rand_op1, &rand_job1, &rand_op2, &rand_job2);
				else {
					macho_op_move(&rand_op1, &rand_job1, &rand_op2, &rand_job2);
					machoflag = TRUE;
				}
			   }
			} 

			if (rand_op1 == rand_op2) {
				current_local_value = schedule_badness(rand_op1);
				swap_job(rand_op1,rand_job1,rand_op2,rand_job2,swapflag);
				new_local_value = schedule_badness(rand_op1);
			} else {
                               	current_local_value = schedule_badness(rand_op1) +
						schedule_badness(rand_op2);
                               	swap_job(rand_op1,rand_job1,rand_op2,rand_job2,swapflag);
                               	new_local_value = schedule_badness(rand_op1) +
						schedule_badness(rand_op2);
			}

			delta = new_local_value - current_local_value;

			flip = random_float(0,1);
			exponent = -(delta)/(K*temperature);
			merit = pow(E,exponent);
			if (merit >= 1.0)
				merit = 0.0; /* don't do unchanging swaps*/

			if (delta < 0) {	/*ACCEPT-WIN*/
				current_value += delta;

				if (TRACE_OUTPUT) {
					if (swapflag) printf("swap "); else printf("move ");
					if (machoflag) printf("macho ");
					printf("WIN %d--%d for operator %d,%d -- delta %d  temp=%f i=%d j=%d\n",
				       rand_job1,rand_job2,rand_op1,rand_op2,delta,temperature,i,j);
			  	}
			} else {
				/*printf("merit=%f flip=%f exponent=%f\n",merit,flip,exponent);*/
				if (merit > flip) { 		/*ACCEPT-LOSS*/
					current_value += delta;
					if (TRACE_OUTPUT) {
					     if (swapflag) printf("swap "); else printf("move ");
					     if (machoflag) printf("macho ");
					     printf("LOSS %d--%d in operators %d,%d -- delta %d  merit=%f flip=%f i=%d j=%d\n",
					 	rand_job1,rand_job2,rand_op1,rand_op2,delta,merit,flip,i,j);
			 		}
				} else { 				/* REJECT */
					swap_job(rand_op2,rand_job2,rand_op1,rand_job1,swapflag);
				}
			}
		}
		if ((current_value - start_value) < 0){	/* rerun at this temp */
			temperature /= COOLING_FRACTION;
			if (TRACE_OUTPUT) printf("rerun at temperature %f\n",temperature);
		}
	}
}

/*	Select a random operator.		*/

int random_operator()
{
	int op;				/* operator to return */

	do {
                op = random_int(0,MAXOPERATIVES-1);
        } while ((operative[op].id == NULL_ID) || (operative[op].jobcount==0));

	return(op);
}

/* Here perform a single operater swap of permuting two jobs on the same operator */

single_op_swap(rand_op1, rand_job1, rand_op2, rand_job2)
int *rand_op1, *rand_job1;	/* return operator and job1 to swap */
int *rand_op2, *rand_job2;      /* return operator and job2 to swap */
{
	/* pick random operative */
	*rand_op2 = *rand_op1 = random_operator();

	/* pick random jobs */
	*rand_job1 = random_int(0,operative[*rand_op1].jobcount-1);
	*rand_job2 = random_int(0,operative[*rand_op2].jobcount-1);
}

/* Here perform a double operater swap of permuting jobs with operators */

double_op_swap(rand_op1, rand_job1, rand_op2, rand_job2)
int *rand_op1, *rand_job1;      /* return operator and job1 to swap */
int *rand_op2, *rand_job2;      /* return operator and job2 to swap */
{
        /* pick random operatives */
	*rand_op1 = random_operator();
	do {
		*rand_op2 = random_operator();
        } while (*rand_op1 == *rand_op2);

        /* pick random jobs */
        *rand_job1 = random_int(0,operative[*rand_op1].jobcount-1);
        *rand_job2 = random_int(0,operative[*rand_op2].jobcount-1);
}

/* Here we perform a move of a job from one operator to another.
   In the current version, we pull a random job from the
   first operator and put it randomly in the second.
*/

double_op_move(rand_op1, rand_job1, rand_op2, rand_job2)
int *rand_op1, *rand_job1;      /* return operator and job1 to move */
int *rand_op2, *rand_job2;      /* return operator and position job2 */
{
        /* pick random operatives */
        *rand_op1 = random_operator();
        do {
                *rand_op2 = random_operator();
        } while ((operative[*rand_op2].id == NULL_ID) ||
		 (overloaded(*rand_op2)) || (*rand_op1 == *rand_op2));

        /* pick random jobs */
        *rand_job1 = random_int(0,operative[*rand_op1].jobcount-1);
        *rand_job2 = random_int(0,operative[*rand_op2].jobcount);
}

/* 	Here we perform a "macho" move of a job from one operator to another.
	Specifically, for a random or the worst operator, we identify the worst
	job in terms of badness, then move it to a random nearby operator.
	We insert it at the place which results in the schedule of least
	badness.
*/

macho_op_move(rand_op1, rand_job1, rand_op2, rand_job2)
int *rand_op1, *rand_job1;      /* return operator and job1 to move */
int *rand_op2, *rand_job2;      /* return operator and position job2 */
{
	int i;			/* counter */
	int index;		/* index of best choice to date */
	int cost;		/* cost of best choice to date */
	int jobindex;		/* absolute job index */
	double random_float();
	
        /* pick either a random operative or the worst operative*/
	if (random_float(0,1) > 0.5)
		*rand_op1 = random_operator();
	else {
		index = 0;
		cost = -LARGE_INT;
		for (i=0; i<MAXOPERATIVES; i++)
			if ((operative[i].id != NULL_ID) && (schedule_badness(i) > cost)) {
				cost = schedule_badness(i);
				index = i;
			}
		*rand_op1 = index;
	}

	/* pick worst job for operative */
	index = 0;
	cost = -LARGE_INT;
	for (i=0; i<operative[*rand_op1].jobcount; i++)
		if (cost <= job_badness(*rand_op1,i)) {
			index = i;
			cost = job_badness(*rand_op1,i);
		}
	*rand_job1 = index;
	jobindex = operative[*rand_op1].jobs[*rand_job1];
	
	/* pick random operator from those closest to job */
	do {
		*rand_op2 = worklist[jobindex].local[random_int(0,NCLOSE-1)];
	} while ((operative[*rand_op2].id == NULL_ID) ||
		(overloaded(*rand_op2)) || (*rand_op1 == *rand_op2));

	/* pick best place to insert job */
	index = operative[*rand_op2].jobcount; /* bias toward end */
	cost = LARGE_INT;
	for (i=0; i<=operative[*rand_op2].jobcount; i++) {
		copy_job(*rand_op1,*rand_job1,*rand_op2,i);
		if (cost >= schedule_badness(*rand_op2)) {
			index = i;
			cost = schedule_badness(*rand_op2);
		}
		delete_job(*rand_op2,i);
	}
	*rand_job2 = index;

/*printf("macho_op_move: rand_op1=%d rand_job1=%d rand_op2=%d rand_job2=%d\n",
	*rand_op1, *rand_job1, *rand_op2, *rand_job2);*/
}

/*	How bad is a given job in a given schedule?		*/

int job_badness(op,i)
int op;				/* operator */
int i;				/* index of job opt badness */
{
	return (TRAVEL_WEIGHT * operative[op].schedule[i].traveltime +
		operative[op].schedule[i].trouble);
}


/*	Swap two jobs from one operative, recomputing the resulting schedule.
	The indicies of the jobs are the current number in the job list for
	the given operator.

	If swapflag in false, the move is a shift of a job from op1 to
	op2, placing it in the job2_nd position of the list.
*/

swap_job(op1, job1, op2, job2, swapflag)
int op1, op2;			/* operator in question */
int job1, job2;			/* indicies of jobs on operators current list */
bool swapflag;			/* TRUE if it is a swap, not a move */
{
	int i;			/* counter */

	if (swapflag == TRUE) {
		swap_ints(&(operative[op1].jobs[job1]), &(operative[op2].jobs[job2]));
		if (op1 == op2)
			fill_in_schedule_for_operative(op1);
		else {
			fill_in_schedule_for_operative(op1);
			fill_in_schedule_for_operative(op2);
		}
	} else if (op1 != op2) {
		/*printf("BEFORE:"); print_op_list(op1); print_op_list(op2);*/
		copy_job(op1, job1, op2, job2);
		delete_job(op1, job1);
		/*printf("AFTER:"); print_op_list(op1); print_op_list(op2);*/
	}
}

/*	Copy job1 from op1 to position job2 for op2.  Do not delete the original copy!	*/

copy_job(op1, job1, op2, job2)
int op1, op2;                   /* operator in question */
int job1, job2;                 /* indicies of jobs on operators current list */
{
	int i;			/* counter */

	for (i=operative[op2].jobcount-1; i>=job2; i--)
		operative[op2].jobs[i+1] = operative[op2].jobs[i];
	operative[op2].jobs[job2] = operative[op1].jobs[job1];
	operative[op2].jobcount++;

	fill_in_schedule_for_operative(op2);
}

/*	Delete the job from the operators list and reschedule.	*/

delete_job(op, job)
int op;				/* operator in question */
int job;			/* indicies of jobs on operators current list */
{
	int i;			/* counter */

	for (i=job+1; i<=operative[op].jobcount; i++)
		operative[op].jobs[i-1] = operative[op].jobs[i];
	operative[op].jobcount--;
	fill_in_schedule_for_operative(op);
}


/*	Print this list of jobs for this operator 	*/

print_op_list(op)
int op;				/* operator */
{
	int i;

	printf("op=%d: ",op);
	for (i=0; i<operative[op].jobcount; i++)
		printf(" %d",operative[op].jobs[i]);
	printf("\n");
}


/****************************************************************************/
/*	UTILITY FUNCTIONS FOR SIMULATED ANNEALING			    */
/****************************************************************************/

/*	These functions generate random numbers in the designated ranges */

random_int(low,high)
int low,high;				/*lower/upper bounds on numb*/
{
	int rand();
	int i,j,r;			/*random number*/

	i = RAND_MAX / (high-low+1);
	i *= (high-low+1);
	while ((j = rand()) >=i) continue;
	r = (j % (high-low+1)) + low;
	if ((r < low) || (r > high))
		printf("ERROR: random integer %d out of range [%d,%d]\n",
			r,low,high);

	/*printf("low=%d high=%d r=%d\n",low,high,r);*/
	return(r);
}

double random_float(low,high)
int low,high;           		/*lower/upper bounds on numb*/
{
	double i,j;			/* avoid arithmetic trouble */
	double r;			/*random number*/
	int rand();

	i = RAND_MAX / (high-low);
        i *= (high-low);
        while ((j = rand()) >=i) continue;
        r = (j/i) * (high-low) + low;

        if ((r < low) || (r > high))
                printf("ERROR: random real %f out of range [%d,%d]\n",
                        r,low,high);

        return(r);
}

