/* objective.cc
 *
 * Simple examples of how to wrap--or otherwise connect--an
 *   objective function so it can work with the DirectSearch
 *   class.
 *
 * Anne Shepherd, 8/2000 (Thanks to Chris Siefert for the EvalF
 *   function).  Revised, pls, 1/2001
 * Virginia Torczon.  Revised 8/2001.
 */

#include "objective.h" 

#define ERROR HUGE_VAL

// two handy little functions that make
// use of the C++ conditional operator ? :

inline double max(double i, double j) {
  return (i > j ? i : j);
} // end max
inline int max(int i, int j) {
  return (i > j ? i : j);
} // end max

// forward declaration
double EvalF(const Vector<double> &point, bool &success, void * VFname);


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
//  fcn() is the default function name for
//    DirectSearch.  The three "fcn" functions below are example of how
//    to call different functions in different ways.
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

/* Here, we use fcn to wrap the user's 
 * objective function.  
 */
void fcn(long vars, Vector<double> &x, double & f, bool & flag, void* nothing)
{
    
  // sanity check 
  if (vars < 1) {
      cerr << "\nError: Dimension cannot be less than 1!!\n";
      exit(1);
  }
  // Here's that actual function call:  change this to reflect what
  // you want to do.
  f = mySin(x);
  // This flag can be used to check for success. Here we don't need it.
  flag = true;
  return;
}

void fcn_2(long vars, Vector<double> &x, double & f,
           bool & flag, void* nothing)
{
     
  if (vars < 1) {
      cerr << "\nError: Dimension cannot be less than 1!!\n";
      exit(1);
  }
  f = myCos(x);
  flag = true;
  return;
}

/**
 * This uses the EvalF() function to call an executable.  It gets its
 * "f" value from the stdout of the objective function.  See user's manual.
 * void * VFname is used to send in the name of the executable.
 */

void fcn_eval(long vars, Vector<double> &x, double & f, bool & flag,
              void * VFname)
{
    
     Vector<double> vec(x);
     f = EvalF(vec, flag, VFname);
}



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
//  The first group of functions are here as simple examples to show how
//    the interface works.
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

// This is a wrapper for the function eval_parab(double x)

void parab_fcn(long vars, Vector<double> & x, double & f,
               bool & flag, void * nothing)
{
  if (vars != 2) {
      cerr << "\nError: Dimension must be 2 !!\n";
      exit(1);
  }
  f = eval_parab(x);
  flag = true;
  return;
}

//  We'll minimize the sine of a number--here vec1
//  is a Vector of one member. (Or even if it isn't, we just
//  take the sine of the first element.)

double mySin(Vector<double> & vec1)
{
    double scalar = vec1.begin()[0];
    return sin(scalar);
}
    

// and here we find the cosine
 
double myCos(Vector<double> & vec2)
{
    double scalar = vec2.begin()[0];
    return cos(scalar);
}


// here, it's an elliptical paraboloid--min at
// the origin. n must equal 2. The wrapper is above.
 
double eval_parab(Vector<double> & parab_vec)
{
    double my_x = parab_vec[0];
    double my_y = parab_vec[1];
    
    return ( (my_x * my_x) + ( 3 * (my_y * my_y) ) );
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
//  Now, a selection of test functions, most of which were
//  coded by Virginia Torczon.
//
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McKinnon @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
// source:  K.I.M. McKinnon, ``Convergence of the Nelder-Mead Simplex
//     Method to a Nonstationary Point,'' Siam J. Optim., Vol 9, No. 1,
//     pp. 148-158. (1998)
//
//   vars must equal 2
//
//   Note that this function is not wrapped--I just wrote it to
//   take the parameters required by the searches.    

void mcKinnon(long vars, Vector<double> & x, double & f,
              bool & flag, void * nothing)
{
  
    double my_x = x[0];
    double my_y = x[1];
    if (my_x <= 0) {
        f = (360 * (my_x * my_x)) + my_y + (my_y * my_y);
    }
    else {
        f = (6 * (my_x * my_x)) + my_y + (my_y * my_y);
    }
    flag = true;
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Avriel 1 @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//    Computes the  test problem (9.B)(i) given in the Exercises for
//    the chapter on Multidimensional Unconstrained Optimization
//    without Derivatives
//    Source:  Mordecai Avriel
//             Nonlinear Programming:  Analysis and Methods.
//
//                            2           2                2
//             f(x) = (x  - 3)  + (x  - 2)  + (x  + x  - 4)
//                      1           2           1    2
//    starting with
//             x  = (0,8),   x  = (0,9),   x  = (1,9)
//              0             1             2
//
//             x  = (8/3,5/3)    with   f(x ) = 1/3
//              *                          *
void avriel1(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

  double f1, f2, f3, x1, x2;
  x1 = x[0];
  x2 = x[1];

  f1 = (x1 - 3.0);
  f2 = (x2 - 2.0);
  f3 = (x1 + x2 - 4.0);

  f = f1*f1 + f2*f2 + f3*f3;
      
  flag=true;
}  // end avriel1

void ComputeGradient1(long vars, Vector<double> & x,
		     Vector<double> & gradient, bool & defined) {

  //  Computes the gradient for test problem (9.B)(i).
  //  Source:  Mordecai Avriel
  //           Nonlinear Programming:  Analysis and Methods.

  //  VARIABLES
  double x1, x2;

  defined = true;

  x1 = x(1);
  x2 = x(2);
  gradient(1) = 4.0*x1 - 14.0 + 2.0*x2;
  gradient(2) = 2.0*x1 - 12.0 + 4.0*x2;

} // end ComputeGradient1 (for Avriel 1)


//@@@@@@@@@@@@@@@@@@@@@@ Ill-Conditioned Avriel 1 @@@@@@@@@@@@@@@@@@@@@@@
//    Computes the  test problem (9.B)(i) given in the Exercises for
//    the chapter on Multidimensional Unconstrained Optimization
//    without Derivatives
//    Source:  Mordecai Avriel
//             Nonlinear Programming:  Analysis and Methods.
//
//    We first observe that the test problem (9.B)(i) can be written in
//    its quadratic form as
//                              T      T
//                  f(x) = 1/2 x Ax + b x + c
//    with
//         A = [4 2    b = [-14    and c = 29.
//              2 4],       -12],
//
//    We then play with the Hessian by replacing A with:
//
//         A = [4 2          where we allow eps -> 0.
//              2 1 + eps],
//
//    The smaller, eps, the more ill-conditioned the resulting
//    Hessian of f(x) since the columns become more and more
//    linearly dependent.  (Obviously, eps = 3 restores us to the
//    problem defined in Avriel.)
//
//    Start the search with
//             x  = (0,8)
//              0
//
//             x  = ((1/2)(7 - 5/eps), 5/eps)
//              *
//   f(x ) = 2*x1*x1 + 2*x1*x2 + (1+eps)*x2*x2 - 14*x1 - 12*x2 + 29
void ill_conditioned_avriel1(long vars, Vector<double> &x,
			     double & f, bool & flag, void* nothing){
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }
  //  VARIABLES
  double x1, x2;

  // Choosing smaller values of eps makes the Hessian
  // ever more ill-conditioned.
  double eps = 0.0001;
  double a11, a12, a21, a22;
  double b1, b2;
  double c;

  a11 = 4.0;
  a12 = 2.0;
  a21 = 2.0;
  a22 = 1.0 + eps;

  b1 = -14.0;
  b2 = -12.0;

  c = 29.0;

  x1 = x[0];
  x2 = x[1];

  f = 0.5*(a11*x1*x1 + a21*x1*x2 + a12*x1*x2 + a22*x2*x2) + b1*x1 + b2*x2 + c;
  
  flag=true;
} // end ill_conditioned_avriel1

void ComputeGradient_ill_conditioned_1(long vars, Vector<double> & x,
				       Vector<double> & gradient, bool & defined) {
  //  Computes the gradient for the redefined test problem (9.B)(i) from
  //  Source:  Mordecai Avriel
  //           Nonlinear Programming:  Analysis and Methods.

  //  VARIABLES
  // Choosing smaller values of eps makes the Hessian
  // ever more ill-conditioned.
  double eps = 0.0001;
  double a11, a12, a21, a22;
  double b1, b2;
  double x1, x2;

  a11 = 4.0;
  a12 = 2.0;
  a21 = 2.0;
  a22 = 1.0 + eps;

  b1 = -14.0;
  b2 = -12.0;

  defined = true;

  x1 = x(1);
  x2 = x(2);
  gradient(1) = a11*x1 + a12*x2 + b1;
  gradient(2) = a21*x1 + a22*x2 + b2;

} // end ComputeGradient_ill_conditioned_1


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Avriel 2 @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//    Computes the test problem (9.B)(ii) given in the Exercises for
//    the chapter on Multidimensional Unconstrained Optimization
//    without Derivatives
//    Source:  Mordecai Avriel
//             Nonlinear Programming:  Analysis and Methods.
//
//                                  2                         2
//             f(x) = (-6 - x  - x )  + (2 - 3x  - 3x  - x x )
//                           1    2            1     2    1 2
//    starting with
//             
//             x  = (-4,6),   x  = (-4,7),   x  = (-3,6)
//              0              1              2
//
//             x  = (0,0),   f(x ) = 40
//              *               *
void avriel2(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

    double f1, f2, x1, x2;

    x1 = x[0];
    x2 = x[1];

    f1 = (-6.0 - x1 - x2);
    f2 = (2.0 - 3.0*x1 - 3.0*x2 - x1*x2);

    f = f1*f1 + f2*f2;

    flag=true;
}  // end avriel2


void ComputeGradient2(long vars, Vector<double> & x,
		     Vector<double> & gradient, bool & defined) {

  //  Computes the gradient for test problem (9.B)(ii).
  //  Source:  Mordecai Avriel
  //           Nonlinear Programming:  Analysis and Methods.

  //  VARIABLES
  double x1, x2;

  defined = true;

  x1 = x(1);
  x2 = x(2);
  gradient(1) = 16.0*x2 + 12.0*x1*x2 +  6.0*x2*x2 + 2.0*x1*x2*x2 + 20.0*x1;
  gradient(2) = 16.0*x1 +  6.0*x1*x1 + 12.0*x1*x2 + 2.0*x1*x1*x2 + 20.0*x2;

} // end ComputeGradient2 (for Avriel 2)



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Beale @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the Beale function:
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (5) Beale function
//
//     (a) n = 2, m = 3
//                              i
//     (b) f (x) = y  - x (1 - x ), where
//          i       i    1      2
//
//         y  = 1.5, y  = 2.25, y  = 2.625
//          1         2          3
//
//     (c) x  = (1,1)
//          0
//
//     (d) f(x ) = 0 at x  = (3, 0.5)
//            *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Beale @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The original version of this test function is:
//
//   E. M. L. Beale
//   On an iterative method of finding a local minimum of a function
//   of more than one variable.
//   Technical Report Number 25,
//   Statistical Techniques Research Group,
//   Princeton University, Princeton, New Jersey.
//   1958.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Beale @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void beale(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

  f=0.0;
  double x1, x2, f1, f2, f3;
  x1 = x[0];
  x2 = x[1];

  f1 = 1.5   - x1 * (1.0 - x2);
  f2 = 2.25  - x1 * (1.0 - x2*x2);
  f3 = 2.625 - x1 * (1.0 - x2*x2*x2);

  f = f1*f1 + f2*f2 + f3*f3;

  flag=true;
}  // end beale


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Box @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the Box three-dimensional function
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (12) Box three-dimensional function
//
//     (a) n = 3, m >= n, variable
//
//     (b) f (x) = exp[-t x ] - exp[-t x ] - x (exp[-t ] - exp[-10t ])
//          i            i 1          i 2     3       i            i
//
//         where t  = (0.1)i
//                i
//
//     (c) x  = (0, 10, 20)
//          0
//
//     (d) f(x ) = 0 at x  = (1, 10, 1), (10, 1, -1)
//            *          *
//          and wherever (x  = x  and x  = 0)
//                         1    2      3
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Box @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   M. J. Box.
//   A comparison of several current optimization methods,
//   and the use of transformations in constrained problems.
//   The Computer Journal,
//   Vol. 9, 1966, pp. 67-77.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Box @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void box(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 3) {
        cerr << "\nThis function is defined for only three variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

  int m = vars;
  double x1, x2, x3, ti, temp1, temp2, temp3, temp4, fi;
  x1 = x[0];
  x2 = x[1];
  x3 = x[2];

  f=0.0;
  for (int i = 1; i <= m; i++) {
    ti = 0.1 * i;
    temp1 = ti*x1;
    temp1 = exp(-temp1);
    temp2 = ti*x2;
    temp2 = exp(-temp2);
    temp3 = exp(-ti);
    temp4 = 10.0*ti;
    temp4 = exp(-temp4);
    fi = temp1 - temp2 - x3*(temp3 - temp4);
    f = f + fi*fi;
  } // end for

  flag=true;
}  // end box

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Brown @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the Brown badly scaled function.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (4) Brown badly scaled function
//
//     (a) n = 2, m = 3
//                        6
//     (b) f (x) = x  - 10
//          1       1   
//                         -6
//         f (x) = x - 2*10
//          2       2
//
//         f  (x) = x x - 2
//          3        1 2
//
//     (c) x  = (1,1)
//          0
//                              6     -6
//     (d) f(x ) = 0 at x  = (10, 2*10  )
//            *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Brown @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is
//   unpublished.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Brown @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void brown(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

  f=0.0;
  double x1, x2, f1, f2, f3;
  x1 = x[0];
  x2 = x[1];

  f1 = x1 - 1000000.0;
  f2 = x2 - 2e-06;
  f3 = x1 * x2 - 2.0;

  f = f1*f1 + f2*f2 + f3*f3;

  flag=true;
}  // end brown


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Helical Valley @@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the Helical Valley function.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (4)  Helical Valley Function
//
//          (a) n = 3
//
//          (b) f (x) = 10(x  - 10 Theta(x ,x ))
//               1          3             1  2
//                           2    2 1/2
//              f (x) = 10((x  + x )    - 1)
//               2           1    2
//
//              f (x) = x
//               3       3
//
//                      3        2
//              f(x) = sum  f (x)
//                     j=1   j
//
//          (c) x  = (-1,0,0)
//               0
//
//          (d) f(x ) = 0 at x  = (1,0,0)
//                 *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Helical Valley @@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   R. Fletcher and Michael J. D. Powell.
//   A rapidly convergent descent method for minimization.
//   The Computer Journal,
//   Vol. 6, No. 2, 1963, pp. 163-168.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Helical Valley @@@@@@@@@@@@@@@@@@@@@@@@
void helical(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if (vars != 3) {
        cerr << "\nThis function is defined for only three variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

    // This is an incredibly sleazy way to initialize pi.
    double pi = 3.14159265358979323846264;
    double factor = 1.0/(2.0*pi);
    double x1, x2, x3, temp, f1, f2, f3;
    f=0.0;
    x1 = x[0];
    x2 = x[1];
    x3 = x[2];

    temp = x1*x1 + x2*x2;

    f1 = 10.0*( x3 - 10.0*(factor*atan(x2/x1)) );
    f2 = 10.0*( sqrt(temp) - 1.0 );
    f3 = x3;

    if (x1 < 0.0) f1 = f1 - 50.0;
      
    f = f1*f1 + f2*f2 + f3*f3;

    flag=true;
}  // end helical

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Penalty I @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes Penalty function I
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (23) Penalty function I
//
//     (a) n variable, m = n + 1
//
//                  1/2
//     (b) f (x) = a   (x  - 1), 1 <= i <= n
//          i            i
//                     n   2
//         f   (x) = (sum x ) - 1/4
//          n+1       j=1  i
//
//                     -5
//         where a = 10
//
//     (c) x  = (epsilon ) where epsilon  = j
//          0           j               j
//
//                             -5
//     (d) f(x ) = 2.24997...10   if n = 4
//            *                -5
//         f(x ) = 7.08765...10   if n = 10
//            *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Penalty I @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   P. E. Gill, W. Murray, and R. A. Pitfield.
//   The implementation of two revised quasi-Newton algorithms for
//   unconstrained optimization.
//   Report NAC 11, National Physics Laboratory.
//   April 1972, pp. 82-83.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Penalty I @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void penalty1(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if ((vars != 4) && (vars != 10)) {
        cerr << "\nThe known minimizer is given only for four and ten variables. ";
        cerr << "\nExiting with value 1.\n";
        exit(1);
    }

    double a = 1.0e-5;
    double sum = 0.0;
    double xJ, fJ;

    f = 0.0;
    for (int j = 0; j < vars; j++) {
      xJ = x[j];
      fJ = xJ - 1.0;

      f = f + a * fJ*fJ;

      sum = sum + xJ*xJ;
    } // end for

    sum = sum - 0.25;
    f = f + sum*sum;

    flag=true;
}  // end penalty1

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Powell @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the EXTENDED Powell singular function any n that
//     is a multiple of 4.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (22) Extended Powell singular function
//
//          (a) n variable but a multiple of 4
//
//          (b) f    (x) = (x     + 10x    )
//               4j-3        4j-3      4j-2
//                          1/2
//              f    (x) = 5  (x     - x  )
//               4j-2           4j-1    4j
//                                         2
//              f    (x) = (x     - 2x    )
//               4j-1        4j-2     4j-1
//                         1/2             2
//              f  (x) = 10   (x     - x  )
//               4j             4j-3    4j
//                     n/4          2           2           2         2
//              f(x) = sum (f    (x)  + f    (x)  + f    (x)  + f  (x) ).
//                     j=1   4j-3        4j-2        4j-1        4j
//
//          (c) x  = (epsilon ) where
//               0           j
//
//          epsilon    = 3, epsilon    = -1, epsilon    = 0, epsilon  = 1
//                 4j-3            4j-2             4j-1            4j
//
//          (d) f(x ) = 0 at x  = (0,...,0)
//                 *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Powell @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original (four dimensional only) version of
//   this test function is:
//
//   M. J. D. Powell.
//   An iterative method for finding stationary values of a
//   function of several variables.
//   The Computer Journal,
//   Vol. 5, 1962, pp. 147-151.
//
//   The source of the extension to higher dimensions is:
//
//   E. Spedicato.
//   Computational experience with quasi-Newton algorithms for
//   minimization problems of moderately large size.
//   Report CISE-N-175,
//   Segrate (Milano).
//   1975.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Powell @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void powell(long vars, Vector<double> &x,
              double & f, bool & flag, void* nothing)
{
    if ((vars < 4) || ((vars%4) != 0)) {
        cerr << "\nThis function requires the number of variables to be a positive multiple of 4. ";
        cerr << "\nExiting with value 1.\n";
        exit(1);
    }


    double f1, f2, f3, f4, x1, x2, x3, x4;

    f = 0.0;
    int top = vars/4;
    for (int j = 0; j < top; j++) {
      x1 = x[4*j];
      x2 = x[4*j+1];
      x3 = x[4*j+2];
      x4 = x[4*j+3];

      f1 = x1 + 10.0*x2;
      f2 = x3 - x4;
      f3 = x2 - 2.0*x3;
      f3 = f3*f3;
      f4 = x1 - x4;
      f4 = f4*f4;

      f = f + f1*f1 + 5.0*f2*f2 + f3*f3 + 10.0*f4*f4;
    } // end for
    flag = true;
} // end powell

//@@@@@@@@@@@@@@@@@@@@@ Powell badly scaled function @@@@@@@@@@@@@@@@@@@@
//     Computes the Powell badly scaled function.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (3) Powell badly scaled function
//
//          (a) n = 2, m = 2
//
//                        4
//          (b) f (x) = 10 * x1 * x2 - 1
//               1
//
//              f (x) = exp[-x1] + exp[-x2] - 1.0001
//               2
//
//          (c) x  = (0,1)
//               0
//                                              -5
//          (d) f(x ) = 0 at x  = (1.098114... 10,  9.106522...)
//                 *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Powell @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original test function is:
//
//   M. J. D. Powell.
//   A hybrid method for nonlinear equations.
//   Numerical Methods for Nonlinear Equations,
//   P. Rabinowitz, editor.
//   Gordon & Breach, New York.
//   1970, pp. 87-114.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Powell @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void powell_badly_scaled(long vars, Vector<double> &x,
              double & f, bool & flag, void* nothing)
{
    if (vars != 2) {
        cerr << "\nThis function is defined for only two variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

    double f1, f2, x1, x2;

    x1 = x[0];
    x2 = x[1];

    f1 = 10000.0 * x1 * x2 - 1.0;
    f2 = exp(-x1) + exp(-x2) - 1.0001;

    f = f1*f1 + f2*f2;

    flag = true;
} // end powell_badly_scaled

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Rosenbrock @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//    Computes the EXTENDED Rosenbrock function for any even value of n.
//
//    Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//             Testing Unconstrained Optimization Software.
//             ACM Transactions on Mathematical Software,
//             Vol. 7, No. 1, March 1981, Pages 17--41.
//
//    (21) Extended Rosenbrock function
//
//         (a) n variable but EVEN
//                                  2
//         (b) f    (x) = 10(x   - x    )
//              2i-1          2i    2i-1
//
//             f  (x) = 1 - x
//              2i           2i-1
//                                                   2         2
//             f(x) = sum from i = 1 to n/2 (f    (x)  + f  (x) ).
//                                            2i-1        2i
//
//         (c) x  = (epsilon ) where epsilon    = -1.2, epsilon  = 1
//              0           i               2i-1               2i
//
//         (d) f(x ) = 0 at x  = (1,...,1)
//                *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Rosenbrock @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original (two dimensional only) version of
//   this test function is:
//
//   H. H. Rosenbrock.
//   An automatic method for finding the greatest
//   or least value of a function.
//   The Computer Journal,
//   Vol. 3, No. 3, October 1960, pp. 175-184.
//
//   The source of the extension to higher dimensions is:
//
//   E. Spedicato.
//   Computational experience with quasi-Newton algorithms for
//   minimization problems of moderately large size.
//   Report CISE-N-175,
//   Segrate (Milano).
//   1975.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Rosenbrock @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void rosenbrock(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
    if ((vars < 2) || ((vars%2) != 0)) {
        cerr << "\nThis function requires a positive, even number of variables. ";
        cerr << "\nExiting with value 1.\n";
        exit(1);
    }

  f=0.0;
  double x1, x2, f1, f2;
  int top = vars/2;
  for(int i=0; i< top; i++) {
    x1 = x[2*i];
    x2 = x[2*i+1];

    f1 = 10.0 * (x2 - x1*x1);
    f2 = 1.0 - x1;

    f += f1*f1 + f2*f2;
  } // end for

  flag=true;
}  // end rosenbrock


//@@@@@@@@@@@@@@@@@@@@@@@ Strictly Differentiable @@@@@@@@@@@@@@@@@@@@@@
//  Computes a strictly differentiable function.
//
//  Source:  Tamara G. Kolda and Robert Michael Lewis
//
//  A strictly differentiable, two-dimensional function
//  that has a minimum f(x ) = 0 at x  = (0,0).
//                         *          *
//
//  Any critical point should be of the form (a, a), where a is any
//  real number, since that is the line along which the function is
//  nondifferentiable---except at the origin.  However, only for a = 0
//  will f(a,a) = 0.  (Of course, numerically, if a is "close" to
//  machine epsilon, 0 <= f(a,a) <= macheps.)
//
void strict(long vars, Vector<double> &x,
	    double & f, bool & flag, void* nothing) {
   
  if (vars != 2) {
    cerr << "\nThis function is defined for only two variables. ";
    cerr << "Exiting with value 1.\n";
    exit(1);
  }

  f=0.0;
  double c11 = 1;
  double c12 = -c11;
  double c21 = -c11;
  double c22 = -c12;
  double x1, x2, xsqr;
  double temp1, temp2, diff1, diff2;

  x1 = x[0];
  x2 = x[1];
  xsqr = x1*x1 + x2*x2;

  temp1 = x1 - c11;
  temp2 = x2 - c12;
  diff1 = temp1*temp1 + temp2*temp2;

  temp1 = x1 - c21;
  temp2 = x2 - c22;
  diff2 = temp1*temp1 + temp2*temp2;

  f = (1.0 - exp(-1000.0 * xsqr)) * 0.5 * max(diff1,diff2);

  flag=true;
}  // end strict


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Trigonometric @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the trigonometric function.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software.
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (3)  Trigonometric function
//
//          (a) n = any positive integer
//                           n
//          (b) f (x) = n - sum (cos x  + j(1 - cos x )  - sin x )
//               j          k=1       k              j          j
//                           n
//                    = n - sum (cos x ) - nj(1 - cos x ) + n sin x
//                          k=1       k                j           j
//                                            2
//              f(x) = sum from j = 1, n f (x)
//                                        j
//          (c) x  = (1/n,...,1/n)
//               0
//
//          (d) f(x ) = 0
//                 *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Trigonometric @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   E. Spedicato.
//   Computational experience with quasi-Newton algorithms for
//   minimization problems of moderately large size.
//   Report CISE-N-175,
//   Segrate (Milano).
//   1975.
//
//   Be warned, however, that there are other variants with the same name.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Trigonometric @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void trigonometric(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {


  double fJ, sum, xJ;

  sum = cos(x[0]);
  for (int j = 1; j < vars; j++) {
    sum = sum + cos(x[j]);
  } // end for

  f = 0.0;
  for (int j = 0; j < vars; j++) {
    xJ = x[j];

    fJ = vars - sum - vars*j*(1.0 - cos(xJ)) + vars*sin(xJ);

    f = f + fJ*fJ;
  } // end for

  flag=true;
}  // end trigonometric


//@@@@@@@@@@@@@@@@@@@@@@@@@ Variably Dimensioned @@@@@@@@@@@@@@@@@@@@@@@@@
//    Computes the Variably dimensioned function.
//
//    Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//             Testing Unconstrained Optimization Software.
//             ACM Transactions on Mathematical Software,
//             Vol. 7, No. 1, March 1981, Pages 17--41.
//
//    (25) Variably dimensioned function
//
//     (a) n variable, m = n + 2
//
//     (b) f (x) = x  - 1, i = 1,...,n
//          i       i
//
//                    n
//         f   (x) = sum j(x  - 1)
//          n+1      j=1    j
//
//                     n            2
//         f   (x) = (sum j(x  - 1))
//          n+2       j=1    j
//
//     (c) x  = (epsilon ) where epsilon  = 1 - (j/n)
//          0           j               j
//
//     (d) f(x ) = 0 at x  = (1,...,1)
//            *          *
//@@@@@@@@@@@@@@@@@@@@@@@@@ Variably Dimensioned @@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is
//   unpublished.
//@@@@@@@@@@@@@@@@@@@@@@@@@ Variably Dimensioned @@@@@@@@@@@@@@@@@@@@@@@@@
void variable(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {

  double fJ, sum;

  sum = 0.0;
  f = 0.0;

  for (int j = 0; j < vars; j++) {
    fJ = x[j] - 1.0;

    f = f + fJ*fJ;

    sum = sum + j*fJ;
  } // end for

  sum = sum*sum;
  f = f + sum;

  sum = sum*sum;
  f = f + sum;
   
  flag=true;
}  // end variable

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Watson @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//     Computes the Watson function.
//
//     Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//              Testing Unconstrained Optimization Software
//              ACM Transactions on Mathematical Software,
//              Vol. 7, No. 1, March 1981, Pages 17--41.
//
//     (26) Watson function
//
//          (a) 2 <= n <= 31, m = 31
//                       n            k-2     n     k-1 2
//          (b) f (x) = sum (k - 1)x t    - (sum x t   )  - 1
//               j      k=2         k j      k=1  k j
//                                            2
//              where t  = j/29, 1 <= j <= 29
//                     j
//
//              f  (x) = x
//               30       1
//                             2
//              f  (x) = x  - x  - 1
//               31       2    1
//
//                      31        2
//              f(x) = sum (f (x))
//                     j=1   j
//
//          (c) x  = (0,...,0)
//               0
//                               -3
//          (d) f  = 2.28767...10    if n = 6
//               *               -6
//              f  = 1.39976...10    if n = 9
//               *               -12
//              f  = 4.72238...10    if n = 12
//               *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Watson @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   J. S. Kowalik and M. R. Osborne.
//   Methods for Unconstrained Optimization Problems.
//   Elsevier, New York.
//   1968.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Watson @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void watson(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
   
    if ((vars != 6) && (vars != 9) && (vars != 12)) {
        cerr << "\nThe known minimizer is given only for six, nine, and twelve variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

    double fJ, f30, f31, sum1, sum2, t[32];

    //     Compute the function value.
    t[0] = 1.0;
    f = 0.0;
    for (int j = 1; j < 30; j++) {
      //        First compute t  and each of its powers up to n - 1
      //                       j
      //                           k
      //        where t(k) = (j/29).
      t[1] = j / 29.0;
      for(int k = 2; k < vars; k++) {
	t[k] = t[1] * t[k-1];
      } // end for

      //                n            k-2               n     k-1
      //        sum1 = sum (k - 1)x t     and  sum2 = sum x t
      //               k=2         k j                k=1  k j
      //
      sum1 = 0.0;
      sum2 = x[0];
      for (int k = 2; k <= vars; k++) {
	sum1 = sum1 + (k - 1.0)*x[k-1]*t[k-2];
	sum2 = sum2 + x[k-1]*t[k-1];
      } // end for

      fJ = sum1 - sum2*sum2 - 1.0;
      f = f + fJ*fJ;
    } // end for

    f30 = x[0];
    f = f + f30*f30;

    f31 = x[1] - x[0]*x[0] - 1.0;
    f = f + f31*f31; 

    flag=true;
}  // end watson

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Wood @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//  Computes the Wood function.
//
//  Source:  J. J. More, B. S. Garbow and K. E. Hillstrom.
//           Testing Unconstrained Optimization Software.
//           ACM Transactions on Mathematical Software,
//           Vol. 7, No. 1, March 1981, Pages 17--41.
//
//  (5)  Wood function
//
//                   2      2           2       2      2           2
//       f(x) = 100(x  - x )  + (1 - x )  + 90(x  - x )  + (1 - x ) 
//                   1    2           1         3    4           3
//
//                                  2           2
//                   + 10.1((1 - x )  + (1 - x ) ) + 19.8(1 - x )(1 - x )
//                                2           4                2       4
//
//                                                      T
//       which has a minimum f(x ) = 0 at x  = (1,1,1,1).
//                              *          *
//                                                    T
//       Standard starting point is x  = (-3,-1,-3,-1).
//                                   0
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Wood @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//   Historical note:
//
//   The source of the original version of this test function is:
//
//   A. R. Colville.
//   A comparative study of nonlinear programming codes.
//   Report 320-2949.
//   IBM New York Scientific Center.
//   1968.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Wood @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void wood(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {
   
    if (vars != 4) {
        cerr << "\nThis function is defined for only four variables. ";
        cerr << "Exiting with value 1.\n";
        exit(1);
    }

    f=0.0;
    double x1, x2, x3, x4, f1, f2, f3, f4, f5, f6, temp;

    x1 = x[0];
    x2 = x[1];
    x3 = x[2];
    x4 = x[3];

    temp = x1*x1 - x2;
    f1 = 100.0 * temp*temp;
    temp = 1.0 - x1;
    f2 = temp*temp;
    temp = x3*x3 - x4;
    f3 = 90.0 * temp*temp;
    temp = 1.0 - x3;
    f4 = temp*temp;
    temp = 1.0 - x2;
    f5 = 1.0 - x4;
    f5 = 10.10 * (temp*temp + f5*f5);
    f6 = 19.80 * (1.0 - x2) * (1.0 - x4);

    f = f1 + f2 + f3 + f4 + f5 + f6;

    flag=true;
}  // end wood

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ xTx @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//                2
//  Computes ||x||
//                2
//
//  f  = 0 at x  = (0,0)
//   *         *
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ xTx @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void xTx(long vars, Vector<double> &x,
		double & f, bool & flag, void* nothing) {

  f = x * x;
  flag=true;
}  // end xTx


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
//  Ugly utility function to allow functions to be specified from
//    outside this file
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Evalf() @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
// Evalf()
//   written by Chris Siefert, slightly modified by Anne Shepherd
//    
// INPUT: The point at which to evaluate the function; a boolean flag;
//     the file name.
// OUTPUT: The function value.
// EFFECT: This function will do the fork/exec dance, sending the
//     program the parameters on its stdin, then after waiting, it will
//     get a  double from the  program's stdout.  If the program exited
//     with an error, it will return ERROR and set success to false.
//     Otherwise success will be true.
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Evalf() @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

double EvalF(const Vector<double> &point, bool &success, void * VFname) {
  int childpid, pipe1[2], pipe2[2];
  const char* Fname;
  long i;
  long P = point.dim();
  double ret_val;
  if (VFname != NULL) {
     Fname = (const char *) VFname;
  }
  FILE *readpipe;
  FILE *writepipe;
  
  long state;
  
  if ((pipe(pipe1) < 0) || (pipe(pipe2) < 0) ) {
    perror("pipe");
    exit(-1);
  }/*end if*/

  if ((childpid = fork()) < 0) {
    perror("fork");
    exit(-1);
  } else if (childpid > 0) {   /*Parent*/
    close(pipe1[0]); close(pipe2[1]);
    
    /* Write to child on pipe1[1], read from child on pipe2[0]. */
    readpipe=fdopen(pipe2[0],"r");
    writepipe=fdopen(pipe1[1],"w");
    if(writepipe==NULL || readpipe==NULL){
      {perror("fdopen");
      cerr<<"MPEF   : Cannot open pipes... exiting\n";abort();}
    }
    /*Output the vector x*/
    for(i=0;i<P;i++) fprintf(writepipe,"%1.8f ", point[i]); 
    fprintf(writepipe,"\n");
    fflush(writepipe);

    /*Read in the returned value, and WAIT*/
    fscanf(readpipe,"%lf", &ret_val);    
    wait(&state);
    /*Close outstanding pipes*/
    fclose(writepipe);fclose(readpipe);
    
    if(WIFEXITED(state)&&(WEXITSTATUS(state)==0)) {
      /*normal exit - no abort, return code 0*/
      success=true;
      return ret_val;     
    }/*end if*/
    else {
      success=false;
      success = true;
      return ERROR;
    }/*end else*/
    
  } else { /*Child*/
    close(pipe1[1]); close(pipe2[0]);
    /* Read from parent on pipe1[0], write to parent on pipe2[1]. */
    dup2(pipe1[0],0); dup2(pipe2[1],1);
    close(pipe1[0]); close(pipe2[1]);

    if (execlp(Fname, Fname, NULL) < 0) {
      perror("execlp");
      abort();
    }
    return ERROR;
    /* Never returns */
  }/*end else*/

  return ERROR;  /*never returns*/
}/*end EvalF*/