# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <ctime>
# include <cmath>

using namespace std;

double best_nearby(double delta[], double point[], double prevbest,
                   int nvars, double f(double x[], int nvars), int *funevals);
int hooke(int nvars, double startpt[], double endpt[], double rho, double eps,
          int itermax, double f(double x[], int nvars));
double r8_abs(double x);
void timestamp(void);
int main();
double rosenbrock(double x[], int n);

double best_nearby(double delta[], double point[], double prevbest,
                   int nvars, double f ( double x[], int nvars ), int *funevals) {
    double ftmp;
    int i;
    double minf;
    double *z;

    z = new double[nvars];

    minf = prevbest;

    for (i = 0; i < nvars; i++) {
        z[i] = point[i];
    }

    for (i = 0; i < nvars; i++) {
        z[i] = point[i] + delta[i];

        ftmp = f ( z, nvars );
        *funevals = *funevals + 1;

        if (ftmp < minf) {
            minf = ftmp;
        } else {
            delta[i] = - delta[i];
            z[i] = point[i] + delta[i];
            ftmp = f ( z, nvars );
            *funevals = *funevals + 1;

            if (ftmp < minf) {
                minf = ftmp;
            } else {
                z[i] = point[i];
            }
        }
    }

    for (i = 0; i < nvars; i++) {
        point[i] = z[i];
    }

    delete [] z;

    return minf;
}

int hooke(int nvars, double startpt[], double endpt[], double rho, double eps,
          int itermax, double f(double x[], int nvars)) {
    double *delta;
    double fbefore;
    int funevals;
    int i;
    int iters;
    int keep;
    double newf;
    double *newx;
    double steplength;
    double tmp;
    bool verbose = true;
    double *xbefore;

    delta = new double[nvars];
    newx = new double[nvars];
    xbefore = new double[nvars];

    for (i = 0; i < nvars; i++) {
        newx[i] = startpt[i];
    }

    for (i = 0; i < nvars; i++) {
        xbefore[i] = startpt[i];
    }

    for (i = 0; i < nvars; i++) {
        if (startpt[i] == 0.0) {
            delta[i] = rho;
        } else {
            delta[i] = rho * r8_abs (startpt[i]);
        }
    }

    funevals = 0;
    steplength = rho;
    iters = 0;
    fbefore = f ( newx, nvars );
    funevals = funevals + 1;
    newf = fbefore;

    while (iters < itermax && eps < steplength) {
        iters = iters + 1;

        if (verbose) {
            cout << "\n";
            cout << "  FUNEVALS, = " << funevals
                 << "  F(X) = " << fbefore << "\n";

            for ( i = 0; i < nvars; i++ ) {
                cout << "  " << i + 1
                     << "  " << xbefore[i] << "\n";
            }
        }
//
//  Find best new point, one coordinate at a time.
//
        for (i = 0; i < nvars; i++) {
            newx[i] = xbefore[i];
        }

        newf = best_nearby ( delta, newx, fbefore, nvars, f, &funevals );
//
//  If we made some improvements, pursue that direction.
//
        keep = 1;

        while (newf < fbefore && keep == 1) {
            for (i = 0; i < nvars; i++) {
//
//  Arrange the sign of DELTA.
//
                if (newx[i] <= xbefore[i]) {
                    delta[i] = - r8_abs (delta[i]);
                } else {
                    delta[i] = r8_abs (delta[i]);
                }
//
//  Now, move further in this direction.
//
                tmp = xbefore[i];
                xbefore[i] = newx[i];
                newx[i] = newx[i] + newx[i] - tmp;
            }

            fbefore = newf;
            newf = best_nearby (delta, newx, fbefore, nvars, f, &funevals);
//
//  If the further (optimistic) move was bad...
//
            if (fbefore <= newf) {
                break;
            }
//
//  Make sure that the differences between the new and the old points
//  are due to actual displacements; beware of roundoff errors that
//  might cause NEWF < FBEFORE.
//
            keep = 0;

            for (i = 0; i < nvars; i++) {
                if (0.5 * r8_abs(delta[i]) < r8_abs(newx[i] - xbefore[i])) {
                    keep = 1;
                    break;
                }
            }
        }

        if (eps <= steplength && fbefore <= newf) {
            steplength = steplength * rho;
            for (i = 0; i < nvars; i++) {
                delta[i] = delta[i] * rho;
            }
        }

    }

    for (i = 0; i < nvars; i++) {
        endpt[i] = xbefore[i];
    }

    delete [] delta;
    delete [] newx;
    delete [] xbefore;

    return iters;
}

double r8_abs(double x) {
    double value;

    if (0.0 <= x) {
        value = x;
    } else {
        value = -x;
    }
    return value;
}

void timestamp(void) {
# define TIME_SIZE 40

    static char time_buffer[TIME_SIZE];
    const struct tm *tm;
    size_t len;
    time_t now;

    now = time(NULL);
    tm = localtime(&now);

    len = strftime(time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm);

    cout << time_buffer << "\n";

    return;
# undef TIME_SIZE
}

int main() {
    timestamp();

    double *endpt;
    double eps;
    int i;
    int it;
    int itermax;
    int nvars = 2;
    double rho;
    double *startpt;
    double value;

    endpt = new double[nvars];
    startpt = new double[nvars];

    cout << "\n";
    cout << "  HOOKE seeks a minimizer of F(X).\n";
    cout << "  Here we use the Rosenbrock function.\n";
//
//  Starting guess for Rosenbrock.
//
    startpt[0] = 0.0;
    startpt[1] = 0.0;

    cout << "\n";
    cout << "  Initial estimate X =\n";
    cout << "\n";
    for (i = 0; i < nvars; i++) {
        cout << "  " << setw(8) << i + 1
             << "  " << setw(14) << startpt[i] << "\n";
    }

    value = rosenbrock(startpt, nvars);

    cout << "\n";
    cout << "  F(X) = " << value << "\n";
//
//  Call HOOKE.
//
    itermax = 5000;
    rho = 0.5;
    eps = 1.0E-06;

    it = hooke(nvars, startpt, endpt, rho, eps, itermax, &rosenbrock);
//
//  Results.
//
    cout << "\n";
    cout << "  Number of iterations taken = " << it << "\n";
    cout << "\n";
    cout << "  X* = \n";
    cout << "\n";
    for (i = 0; i < nvars; i++) {
        cout << "  " << setw(8) << i + 1
             << "  " << setw(14) << endpt[i] << "\n";
    }

    value = rosenbrock(endpt, nvars);

    cout << "\n";
    cout << "  F(X*) = " << value << "\n";

    delete [] endpt;
    delete [] startpt;

    cout << "  Normal end of execution.\n";

    cout << "\n";
    timestamp();

    return 0;
}

double rosenbrock(double x[], int n) {

    return x[0]*x[0] - 2*sin(x[1]) + 1.85*x[1]*x[1];

}