#define USER_ID \
"/* $Id: user.c,v 7.3 1993/02/04 18:33:47 ingber Exp ingber $ */"

#include "user.h"

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	main - A sample optimization code that calls vfsr
* Parameters: None
* Inputs:
* Outputs:
* Global Variables: USER_OPTIONS [optional]
* Returns: cost_value
* Calls:
*	fopen                 - to open the output file
*	print_time("start")   - prints time
*	initialize_rng() - to initialize the random number generator
*	initialize_parameters - to initialize the
*				parameters and allocate storage space
*	vfsr - the VFSR optimization call to minimize the function
*	print_time - to print the execution time statistics
*	fclose                - close output files
*	exit                  - exit successfully
* Description:
*	This is a sample calling program to optimize using VFSR
* Local Variables:
*	exit_code             - exit code returned by vfsr
*	compile_cnt - number of compile options set by user in Makefile
* Portability:
* Other:
***********************************************************************/
#if HAVE_ANSI
int main(int argc, char **argv)
#else
int main(argc, argv)
int argc;
char **argv;
#endif
{
    int exit_code;
    int compile_cnt;

    /* open the output file */
    ptr_out = fopen("user_out", "w");
    fprintf(ptr_out, "%s\n\n", USER_ID);

    /* print out compile options set by user in Makefile */
    if (argc > 0) {
	fprintf(ptr_out, "CC = %s\n", argv[1]);
	for (compile_cnt = 2; compile_cnt < argc; ++compile_cnt) {
	    fprintf(ptr_out, "\t%s\n", argv[compile_cnt]);
	}
	fprintf(ptr_out, "\n");
    }
#if TIME_CALC
    /* print starting time */
    print_time("start");
#endif

    initialize_rng();

    /* initialize the users parameters, allocating space, etc.
       Note that the default is to have vfsr generate the initial
       cost_parameters that satisfy the user's constraints. */

    initialize_parameters();

    /* optimize the cost_function, returning the results in
       cost_value and cost_parameters */
    cost_value = vfsr(		/* return cost_value */
			 cost_function,	/* (double) (*cost_function)
				  (double *cost_parameters,
                                   double *parameter_lower_bound,
                                   double *parameter_upper_bound,
				   int *cost_flag,
                                   USER_DEFINES *USER_OPTIONS) */
			 randflt,	/* (double) (*random_generator) () on
				   {0,1} this might be drand48 for quick
				   and dirty unix systems */
			 parameter_dimension, parameter_int_real,
			 cost_parameters,	/* array of initial parameters,
					   replaced by final parameters on
					   return of vfsr() */
			 parameter_lower_bound, parameter_upper_bound,
			 cost_tangents, cost_curvature,
			 &exit_code,
			 USER_OPTIONS);

#if TIME_CALC
    /* print ending time */
    print_time("end");
#endif

    /* close all files */
    fclose(ptr_out);

    exit(0);
    /* NOTREACHED */
}

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	initialize_parameters - sample parameter initialization function
* Parameters:
* Inputs:
* Outputs:
* Global Variables:
*	parameter_dimension   - the number of free variables
*	parameter_lower_bound - the minimum parameter value allowed
*	parameter_upper_bound - the maximum parameter value allowed
*	cost_parameters - the parameter initial and final values
*	parameter_int_real - the parameter types, real or integer
*	cost_tangents - parameter derivatives used during reannealing
*                     - can contain initial parameter temperatures
*	cost_curvature - parameter 2nd derivatives
*                     - [0] member can contain initial cost temperature
* Returns: None
* Calls:
*	calloc                - to allocate storage
* Description:
*	This depends on the users cost function to optimize (minimum).
*	The routine allocates storage needed for vfsr. The user should
*	define the number of parameters and their ranges,
*	and make sure the initial parameters are within
*	the minimum and maximum ranges. The array
*	parameter_int_real should be REAL_TYPE (0) for real parameters,
*	and INTEGER_TYPE (1) for integer values
* Local Variables:
*       exit_code             - tells how vfsr exited
* Portability:
* Other:
***********************************************************************/
#if HAVE_ANSI
void initialize_parameters(void)
#else
void initialize_parameters()
#endif
{
    /* the number of parameters for the cost function */
    parameter_dimension = 4;

    /* ALLOCATE STORAGE */
    /* allocate parameter minimum space */
    parameter_lower_bound =
	(double *) calloc(parameter_dimension, sizeof(double));
    /* allocate parameter maximum space */
    parameter_upper_bound =
	(double *) calloc(parameter_dimension, sizeof(double));
    /* allocate parameter initial values; the parameter final values will
       be stored here later */
    cost_parameters =
	(double *) calloc(parameter_dimension, sizeof(double));
    /* allocate the parameter types, real or integer */
    parameter_int_real =
	(int *) calloc(parameter_dimension, sizeof(int));
    /* allocate space for the parameter cost_tangents - used for reannealing */
    cost_tangents =
	(double *) calloc(parameter_dimension, sizeof(double));
    /* allocate space for the parameter cost_curvatures/covariance matrix */
    cost_curvature =
	(double *) calloc(parameter_dimension * parameter_dimension,
			  sizeof(double));
#if ADAPTIVE_OPTIONS
    USER_OPTIONS =
	(USER_DEFINES *) calloc(1, sizeof(USER_DEFINES));
#endif

    /* store some parameter ranges */
    parameter_lower_bound[0] = parameter_lower_bound[1] =
	parameter_lower_bound[2] = parameter_lower_bound[3] = -10000.0;
    parameter_upper_bound[0] = parameter_upper_bound[1] =
	parameter_upper_bound[2] = parameter_upper_bound[3] = 10000.0;

    /* store the initial parameter values */
    cost_parameters[0] = 999.0;
    cost_parameters[1] = -1007.0;
    cost_parameters[2] = 1001.0;
    cost_parameters[3] = -903.0;

    /* store the initial parameter types */
    parameter_int_real[0] = parameter_int_real[1] =
	parameter_int_real[2] = parameter_int_real[3] = REAL_TYPE;

/* These can be used if set DEFINE_OPTIONS = ADAPTIVE_OPTIONS=TRUE
   in the Makefile.  These can be overridden by calls to data file,
   permitting tuning/debugging without requiring recompiling. */

#if ADAPTIVE_OPTIONS
#if OPTIONS_FILE
    ptr_options = fopen("options_file", "r");
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).COST_PRECISION = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).USER_INITIAL_PARAMETERS = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).ACCEPTED_TO_GENERATED_RATIO = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).LIMIT_ACCEPTANCES = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).TEMPERATURE_RATIO_SCALE = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).TEMPERATURE_ANNEAL_SCALE = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).COST_PARAMETER_SCALE = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).TESTING_FREQUENCY_MODULUS = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).MAXIMUM_REANNEAL_INDEX = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).REANNEAL_RESCALE = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).INITIAL_PARAMETER_TEMPERATURE = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).USER_INITIAL_PARAMETERS_TEMPS = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).USER_INITIAL_COST_TEMP = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).NUMBER_COST_SAMPLES = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).MAXIMUM_COST_REPEAT = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%lf", &read_double);
    (*USER_OPTIONS).DELTA_X = read_double;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).INCLUDE_INTEGER_PARAMETERS = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).ACTIVATE_REANNEAL = read_int;
    fscanf(ptr_options, "%s", read_option);
    fscanf(ptr_options, "%d", &read_int);
    (*USER_OPTIONS).LIMIT_INVALID_GENERATED_STATES = read_int;
    fclose(ptr_options);
#else				/* OPTIONS_FILE */
    (*USER_OPTIONS).COST_PRECISION = (double) SMALL_FLOAT;
    (*USER_OPTIONS).USER_INITIAL_PARAMETERS = FALSE;
/*	(*USER_OPTIONS).ACCEPTED_TO_GENERATED_RATIO = 1.0E-6; */
    (*USER_OPTIONS).ACCEPTED_TO_GENERATED_RATIO = 1.0E-4;
/*	(*USER_OPTIONS).LIMIT_ACCEPTANCES = 10000; */
    (*USER_OPTIONS).LIMIT_ACCEPTANCES = 1000;
    (*USER_OPTIONS).TEMPERATURE_RATIO_SCALE = 1.0E-5;
    (*USER_OPTIONS).TEMPERATURE_ANNEAL_SCALE = 100.0;
    (*USER_OPTIONS).COST_PARAMETER_SCALE = 1.0;
    (*USER_OPTIONS).TESTING_FREQUENCY_MODULUS = 100;
    (*USER_OPTIONS).MAXIMUM_REANNEAL_INDEX = 100000000;
    (*USER_OPTIONS).REANNEAL_RESCALE = 10.0;
    (*USER_OPTIONS).INITIAL_PARAMETER_TEMPERATURE = 1.0;
    (*USER_OPTIONS).USER_INITIAL_PARAMETERS_TEMPS = FALSE;
    (*USER_OPTIONS).USER_INITIAL_COST_TEMP = FALSE;
    (*USER_OPTIONS).NUMBER_COST_SAMPLES = 5;
    (*USER_OPTIONS).MAXIMUM_COST_REPEAT = 5;
    (*USER_OPTIONS).DELTA_X = 0.001;
    (*USER_OPTIONS).INCLUDE_INTEGER_PARAMETERS = FALSE;
    (*USER_OPTIONS).ACTIVATE_REANNEAL = TRUE;
    (*USER_OPTIONS).LIMIT_INVALID_GENERATED_STATES = 1000;
#endif				/* OPTIONS_FILE */
#endif				/* ADAPTIVE_OPTIONS */

}

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	double cost_function - a sample cost function
* Parameters:
* Inputs:
*	double *x                     - the array of function parameters
*       double *parameter_lower_bound - array of lower bounds
*       double *parameter_upper_bound - array of upper bounds
* Outputs:
*	int *cost_flag - flag indicating parameters generated were valid
*       double *parameter_lower_bound - array of lower bounds
*       double *parameter_upper_bound - array of upper bounds
* Global Variables:
*       USER_OPTIONS [optional]
* Returns:
*	the user cost value
* Calls:
*	fabs
* Description:
*	This is the users cost function to optimize
*	(find the minimum).
*	cost_flag is set to true if the parameter set
*	does not violates any constraints
*       parameter_lower_bound and parameter_upper_bound may be
*       adaptively changed during the search.
* Local Variables:
*	double ans           - the current results
* Portability:
* Other:
***********************************************************************/
#if HAVE_ANSI
double cost_function(double *x,
		     double *parameter_lower_bound,
		     double *parameter_upper_bound,
		     int *cost_flag,
		     USER_DEFINES * USER_OPTIONS)
#else
double cost_function(x,
		     parameter_lower_bound,
		     parameter_upper_bound,
		     cost_flag,
		     USER_OPTIONS)
double *x;
double *parameter_lower_bound;
double *parameter_upper_bound;
int *cost_flag;
USER_DEFINES *USER_OPTIONS;
#endif
{
/* Objective function from
%A A. Corana
%A M. Marchesi
%A C. Martini
%A S. Ridella
%T Minimizing multimodal functions of continuous variables
with the "simulated annealing" algorithm
%J ACM Trans. Mathl. Software
%V 13
%S 3
%P 272-280
%D 1987

We thank Mark Johnson <mark@microunity.com> for bringing this
function to our attention.

This function contains 1.0E20 local minima.  Visiting each minimum
for a millisecond would take about the present age of the universe
to visit all these minima. */

#define ARRAY_LIMIT	4

    double s[ARRAY_LIMIT], t[ARRAY_LIMIT], d[ARRAY_LIMIT];
    double term, summ, dj, c, deviate;
    int i;
    static long int funevals = 0;

    for (i = 0; i < ARRAY_LIMIT; ++i) {
	s[i] = 0.2;
	t[i] = 0.05;
	d[i] = 10.0;
    }

    d[0] = 1.0;
    d[1] = 1000.0;
    d[2] = 10.0;
    d[3] = 100.0;
    c = 0.15;

    summ = 0.0;
    for (i = 0; i < 4; ++i) {
	dj = floor(fabs(x[i] / s[i]) + 0.49999999);
	if (x[i] < 0.0)
	    dj = 0.0 - dj;
	dj = dj * s[i];
	deviate = fabs(x[i] - dj);
	if (deviate < fabs(t[i])) {
	    term = 0.0;
	    if (dj < 0.0)
		term = t[i];
	    if (dj > 0.0)
		term = (0.0 - t[i]);
	    term = term + dj;
	    term = term * term * c * d[i];
	} else
	    term = d[i] * (x[i] * x[i]);
	summ = summ + term;
    }
    funevals = funevals + 1;

    *cost_flag = TRUE;

#if TIME_CALC
    /* print the time every PRINT_FREQUENCY evaluations */
    if ((PRINT_FREQUENCY > 0) &&
	((funevals % (long int) PRINT_FREQUENCY) == 0)) {
	fprintf(ptr_out, "funevals = %ld  ", funevals);
	print_time("");
    }
#endif
    return (summ);
}

/* ************************************************* */
/* Here is a good random number generator
   and the auxiliary routines to print the execution time */

#define SHUFFLE 256		/* size of random array */
#define MULT 25173
#define MOD ((long int) 65536)
#define INCR 13849
#define FMOD ((double) 65536.0)

static long int seed = 696969;
double random_array[SHUFFLE];	/* random variables */

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	double myrand(void) - returns a random number between 0 and 1
* Parameters: None
* Inputs:
* Outputs:
* Global Variables:
*	double seed 	    - The rng seed
* Returns: None
* Calls: None
* Description:
*	This routine returns the random number generator between 0 and 1
* Local Variables: None
* Portability:
* Other:
***********************************************************************/

#if HAVE_ANSI
double myrand(void)
#else
double myrand()
#endif
/* returns random number in {0,1} */
{
    seed = (MULT * seed + INCR) % MOD;
    return ((double) seed / FMOD);
}

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	double randflt(void) - Shuffles random numbers in random_array[]
* Parameters: None
* Inputs:
* Outputs:
* Global Variables: None
* Returns: None
* Calls:
*	myrand()
* Description:
*	This routine initializes the random number generator
* Local Variables: None
* Portability:
* Other:
***********************************************************************/

#if HAVE_ANSI
double randflt(void)
#else
double randflt()
#endif
/* shuffles random numbers in random_array[SHUFFLE] array */
{
    double rranf;
    int kranf;

    kranf = (int) (myrand() * SHUFFLE) % SHUFFLE;
    rranf = *(random_array + kranf);
    *(random_array + kranf) = myrand();

    return (rranf);
}

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	initialize_rng() - to initialize the random number generator
* Parameters: None
* Inputs:
* Outputs:
* Global Variables: None
* Returns: None
* Calls:
*	myrand()
*	randflt()
* Description:
*	This routine initializes the random number generator
* Local Variables: None
* Portability:
* Other:
***********************************************************************/

#if HAVE_ANSI
void initialize_rng(void)
#else
void initialize_rng()
#endif
{
    int n;
    double x;

    for (n = 0; n < SHUFFLE; ++n)
	random_array[n] = myrand();
    for (n = 0; n < 1000; ++n)	/* warm up random generator */
	x = randflt();
}

#if TIME_CALC
/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	print_time     - print the time routine
* Parameters: None
* Inputs:
*	char *message  - the output message
* Outputs:
* Global Variables: None
* Returns: None
* Calls:
*	getrusage      - get the system resource usage information
*	aux_print_time - to calculate and print the runtime
* Description:
*	This calculates the time and runtime and prints it.
* Local Variables:
*	int who        - name of person (self) to get information on
*	struct rusage usage - holds the resource information
* Portability:
* Other:
***********************************************************************/
#if HAVE_ANSI
int getrusage(int who, struct rusage *usage);

void print_time(char *message)
#else
int getrusage();

void print_time(message)
char *message;

#endif
{
    int who = RUSAGE_SELF;	/* Check our own time */
    struct rusage usage;

    /* get the resource usage information */
    getrusage(who, &usage);

    /* print the usage time in reasonable form */
    aux_print_time(&usage.ru_utime, message);
}

/***********************************************************************
* Author: Lester Ingber, Bruce Rosen (copyright) (c)
* Date	 5 Nov 92
* Procedure:
*	aux_print_time     - auxiliary print the time routine
* Parameters: None
* Inputs:
*	struct timeval *time - the current time
*	char *message      - the output message
* Outputs:
* Global Variables: None
* Returns: None
* Calls: None
* Description:
*	This calculates the time and runtime and prints it.
* Local Variables:
*	static double sx   - previous time in seconds
*	double us, s, m, h - current time in microsec, sec, min, hours
*	double ds, dm, dh  - last time difference in sec, min, hours
* Portability:
* Other:
***********************************************************************/
#if HAVE_ANSI
void aux_print_time(struct timeval *time, char *message)
#else
void aux_print_time(time, message)
struct timeval *time;
char *message;
#endif
{
    static double sx;
    double us, s, m, h;
    double ds, dm, dh;

    /* calculate the new microseconds, seconds, minutes, hours
       and the differences since the last call */
    us = (double) ((int) ((double) SMALL_FLOAT + time->tv_usec)) / 1.E6;
    s = (double) ((int) ((double) SMALL_FLOAT + time->tv_sec)) + us;
    ds = s - sx;
    sx = s;

    h = (int) ((double) SMALL_FLOAT + s / 3600.);
    m = (int) ((double) SMALL_FLOAT + s / 60.) - 60. * h;
    s -= (3600. * h + 60. * m);
    dh = (int) ((double) SMALL_FLOAT + ds / 3600.);
    dm = (int) ((double) SMALL_FLOAT + ds / 60.) - 60. * dh;
    ds -= (3600. * dh + 60. * dm);

    /* print the statistics */
    fprintf(ptr_out,
	    "%s:time: %gh %gm %gs; incr: %gh %gm %gs\n",
	    message, h, m, s, dh, dm, ds);
}

#endif				/* TIME_CALC */