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

extern int cpattrn(int choice, void (*obj)(double *, double *), int n,
                   double x[], double step1, double step2, int maxf,
                   double aux[], int *iaux);
extern int cshh(void (*obj)(double *, double *), int n, double x[],
                int right, double step1[], double step2, int maxf,
                int stddev, double aux[], int *iaux);
extern int cnm(void (*obj)(double *, double *), int n, double x[],
               int right, double sigma, double alpha, double beta,
               double gamma, double step1[], double step2, int maxf,
               int stddev, double aux[], int *iaux);
extern int csmd(void (*obj)(double *, double *), int n, double x[],
                int right, double step1[], double step2, int maxf,
                int stddev, double aux[], int *iaux);

void objc(x, f, n)
     double x[], *f;
     long *n;
{
  /************************************************************************

     This subroutine implements the objective function. Bounds are
     constructed by making the function evaluation very large when
     the x values are out of range.

   ***********************************************************************/
  double temp1, temp2, temp3;
  if (x[0]<6.0 || x[1] < 0.0)
    *f = 10000000.0;
  else {
    temp1 = 1.0 / (27.0 * sqrt(3.0));
    temp2 = x[0] - 3.0;
    temp3 = x[1] * x[1] * x[1];
    *f  = temp1 * (temp2*temp2 - 9.0) * temp3;
  }
}

int mymain_()
{
/***********************************************************************

     This is an example of an interactive program which performs
     a direct search on the objective function "objf". You can do
     a search and replace to substitute your own objective function
     and recompile.
     You can also redirect a file to standard input to avoid user
     interaction:
         MySearch < sample.dat

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

  const int maxn=10;
  int type, i, set, iright, istat=0, istd, right, stddev;
  int n, maxf, iaux = 0;
  double start[maxn], step1[maxn], step2, aux[maxn+2];
  double sigma, alpha, beta, gamma;

  printf("Enter the type of search you would like to perform.\n");
 invalid:
  printf("  1 - Coordinate Search\n");
  printf("  2 - Compass Search\n");
  printf("  3 - NLess Search\n");
  printf("  4 - Hooke and Jeeves\n");
  printf("  5 - Hooke and Jeeves (edited by E. Dolan)\n");
  printf("  6 - Simplex based on Spendley, Hext, and Himsworth\n");
  printf("  7 - Simplex based on Nelder and Meade\n");
  printf("  8 - Sequential Multidirectional Search\n");
  scanf("%d", &type);
  if (type < 1 || type > 8) {
    printf("That was not a valid entry. Choose a number between 1 and 8\n");
    goto invalid;
  }
  printf("Enter the dimension of the problem.\n");
  printf("If you are running the default subroutine objc, then the dimension must be 2.\n");
  scanf("%d", &n);
  if (n > maxn) {
    printf("Increase the constant maxn before proceeding.\n");
    exit(-1);
  }
  printf("Enter the starting values for x.\n");
  for (i=0;i<n;i++)
    scanf("%lf", &start[i]);
  printf( "Enter the maximum number of function calls. (-1 for no max)\n");
  scanf("%d", &maxf);
  /* Input for pattern searches */
  if (type >= 1 && type <= 5) { 
    printf("Enter the initial step length. Enter -1 to get the default of 0.25\n");
    scanf("%lf", &step1[0]);
    printf("Enter the final step length. Enter -1 to get the default of 10E-8\n");
    scanf("%lf", &step2);
    istat = cpattrn(type, objc, n, start, step1[0], step2, maxf, aux, &iaux);
  }
  /* Input for simplex searches */
  else {
    printf("Enter the simplex type: 1-right 2-regular\n");
    scanf("%d", &iright);
    if (iright == 1) 
      right = 1;
    else if (iright==2)
      right = 0;
    else {
      printf("Input was invalid. Right simplex was chosen.\n");
      right = 1;
    }
    printf("Do you want to set the edgelengths or use the default of 2.0? 1-Set 2-Use Default\n");
    scanf("%d", &set);
    if (set == 1) {
      if (right==1) {
	printf("Enter the n edgelengths. (Direct search will calculate length n+1.)\n");
	for (i=0;i<n;i++)
	  scanf("%lf", &step1[i]);
      }
      else {
	printf("Enter one edgelength of the regular simplex.\n");
	scanf("%lf", &step1[0]);
      }
    }
    else if (set==2) {
      step1[0] = -1.0;
    }
    else {
      printf("Input was invalid. Default edgelengths used\n");
      step1[0] = -1.0;
    }
    printf("Stopping criteria on standard deviation or delta?\n");
    printf("1-stddev 2-delta\n");
    scanf("%d", &istd);
    if (istd == 1) 
      stddev = 1;
    else if (istd == 2)
      stddev = 0;
    else {
      printf("Input was invalid. Stopping on stddev.\n");
      stddev = 1;
    }
    printf("Enter the final step length. Enter -1 to get the default of 10E-8\n");
    scanf("%lf", &step2);
    if (type==6)
      istat = csmd(objc, n, start, right, step1, step2, maxf, stddev, aux, &iaux);
    else if (type==7) {
      printf("Enter the shrinking coefficient, sigma.\n");
      printf("(-1 uses a default of 0.5.)\n");
      scanf("%lf", &sigma);
      printf("Enter the reflection coefficient, alpha.\n");
      printf("(-1 uses a default of 1.0.)\n");
      scanf("%lf", &alpha);
      printf("Enter the contraction coefficient, beta.\n");
      printf("(-1 uses a default of 0.5.)\n");
      scanf("%lf", &beta);
      printf("Enter the expansion coefficient, gamma.\n");
      printf("(-1 uses a default of 2.0.)\n");
      scanf("%lf", &gamma);
      istat = cnm(objc, n, start, right, sigma, alpha, beta, gamma, step1, 
		  step2, maxf, stddev, aux, &iaux);
    } 
    else {
      istat = csmd(objc, n, start, right, step1, step2, maxf, stddev, aux,
		   &iaux); 
    }
  }
  if (istat == 1)
  {
    printf("The final delta is: %f \n", aux[0]);
    printf("The minimum function evaluation is: %f \n", aux[1]);
    for (i=2; i<=n+1;i++)
      printf("The corresponding x is: %f \n", aux[i]);
    printf("The number of evaluations was: %d\n", iaux);
  }

  return 0;

}