/*************************************************************************/
/*** DHC.C                                                             ***/
/*** Example code demonstrating the usage of the dynamic hill climbing ***/
/*** algorithm on DeJong's test function Rosenbrock's saddle.	       ***/
/*** Note: compile with 'gcc dhc.c -lm'            		       ***/
/*** Last modified Sep 24, 1995.  (c) Deniz Yuret, deniz@ai.mit.edu    ***/
/*************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* NDIM and THRESHOLD are the only two parameters to experiment with.  */
/* Try increasing NDIM to see the performance of DHC in higher dims.   */
/* Try making the THRESHOLD smaller to get more accurate results, and  */
/*   observe the performance trade-off.                                */
/* If you make THRESHOLD larger, DHC might fail to find the optimum    */
/*   at all, i.e. you may get a qualitatively different behavior.      */
/* Try multiple restarts by giving the program a command line argument.*/
/* After you have played around enough, you should write your own      */
/*   restart method to take advantage of your specific landscape.      */
/* You should take a look at my thesis for details.                    */
/*   (Yuret, MIT MS 1990) ftp://ftp.ai.mit.edu/pub/DHC/dhc-msthesis.ps */
/* Compare your favorite algorithm with the same function/parameters.  */
/*   Send me e-mail if it can beat DHC :)                              */
/* Good luck optimizing.                                               */

#define NDIM 2			/* Number of dimensions */
#define THRESHOLD 1E-4		/* Minimum step size */

#define INIT_SIZE 1.0		/* Initial step size */
int count;			/* Number of function calls */

double func(x)			/* Objective function: */
     double x[];		/* Generalized Rosenbrock's Saddle */
{				/* From DeJong (PhD, U of Mich, 75) */
  int i;
  double sum;

  count++;			/* Increment eval count (no cheating :) */
  for(i=0, sum=0; i<NDIM-1; i++)
    sum += 100*(x[i+1]-x[i]*x[i])*(x[i+1]-x[i]*x[i]) + (x[i]-1)*(x[i]-1);
  return sum;
}

void main(argc,argv)		/* argv[1]: number of restarts */
     int argc;
     char **argv;
{
  void local_optimize();
  double x[NDIM];
  int i,j;
  int try;

  if(argc < 2) try = 1;
  else try = atoi(argv[1]);
  count = 0;
  srand(time(NULL));
  for(j=0; j<NDIM; j++)
    x[j] = -1;
  for(i=0; i<try; i++) {
    printf("\ntry=%d\n\n", i);

    /* Random jump for each restart.                        */
    /* You should do something more intelligent based on    */
    /* the landscape of your problem and previous results   */

    for(j=0; j<NDIM; j++)
      x[j] *= exp((double)rand()/pow(2,31) - 0.5);

    /* local_optimize directly writes the array */
    /* argument x and leaves the answer there   */

    local_optimize(func, x);
  }
}

/* You should really look at my thesis to understand the following code. */

void local_optimize(f,x)
     double (*f)();
     double x[];
{
  double u[NDIM], v[NDIM], xv[NDIM];
  double fx,fxv;
  int vi,vvec;
  double vr;
  int i,iter,maxiter;

  for(i=0; i<NDIM; i++)
    u[i] = v[i] = 0;
  vi = -1; vvec = 1;
  vr = -INIT_SIZE;
  fx = f(x);
  fxv = 1E10;

  printf("%d. %.4f <= ", count, fx);
  for(i=0; i<NDIM; i++) { printf("%.4f ", x[i]); }; printf("\n");

  while(fabs(vr) >= THRESHOLD) {
    maxiter = ((fabs(vr) < 2*THRESHOLD) ? 2*NDIM : 2);
    iter = 0;
    while((fxv >= fx) && (iter < maxiter)) {
      if(iter == 0) { for(i=0; i<NDIM; i++) xv[i] = x[i]; }
      else xv[vi] -= vr;
      if(vvec) vvec = 0;
      vr = -vr;
      if(vr > 0) vi = ((vi+1) % NDIM);
      xv[vi] += vr;
      fxv = f(xv);
      iter++;
    }
    if(fxv >= fx) {
      fxv = 1E10;
      vr /= 2;
    } else {
      fx = fxv; printf("%d. %.4f <= ", count, fx);
      for(i=0; i<NDIM; i++) { x[i] = xv[i]; printf("%.4f ", x[i]); }
      printf("\n");
      if(iter == 0) {
	if(vvec) {
	  for(i=0; i<NDIM; i++) {
	    u[i] += v[i]; v[i] *= 2; xv[i] += v[i];
	  }
	  vr *= 2;
	} else {
	  u[vi] += vr; vr *= 2; xv[vi] += vr;
	}
	fxv = f(xv);
      } else {
	for(i=0; i<NDIM; i++) xv[i] += u[i];
	xv[vi] += vr;
	fxv = f(xv);
	if(fxv >= fx) {
	  for(i=0; i<NDIM; i++) { u[i] = 0; xv[i] = x[i]; }
	  u[vi] = vr; vr *= 2;
	  xv[vi] += vr; fxv = f(xv);
	} else {
	  for(i=0; i<NDIM; i++) x[i] = xv[i]; fx = fxv;
	  u[vi] += vr;
	  for(i=0; i<NDIM; i++) v[i] = 2*u[i]; vvec = 1;
	  for(i=0; i<NDIM; i++) xv[i] += v[i]; fxv = f(xv);
	  for(vr=0,i=0; i<NDIM; i++) vr += v[i]*v[i];
	  vr = sqrt(vr);
	}
      }
    }
  }
}