#include <iostream>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>

double my_f(const gsl_vector *x, void *params);


double my_f(const gsl_vector *x, void *cur_model)
{
    double objf;

    objf = 1.0;

    return(objf);
}

int main(void)
{
    int n = 50;
    //simplex optimization
    /* Set up the starting point */

    const gsl_multimin_fminimizer_type *T =
        gsl_multimin_fminimizer_nmsimplex2;
    gsl_multimin_fminimizer *s=NULL;
    gsl_vector *x, *ss;
    gsl_multimin_function minex_func;

    // initialize x, ss
    x = gsl_vector_alloc(n);
    ss = gsl_vector_alloc (n);

    for (int i=0; i<int(n); i++) {
        gsl_vector_set(x, i, i+0.1);
        gsl_vector_set(ss, i, 0.2);
    }

    std::cout << std::endl;
    std::cout << "initial x vector is" << std::endl;
    for (int i=0; i<int(n); i++) {
        std::cout << gsl_vector_get(x, i) << "  ";
    }
    std::cout << std::endl;
    std::cout << "initial step size vector is" << std::endl;
    for (int i=0; i<int(n); i++) {
        std::cout << gsl_vector_get(ss, i) << "  ";
    }
    std::cout << std::endl;

    double par[1] = {1.0};

    /* Initialize the method */
    minex_func.n = n;
    minex_func.f = &my_f;
    minex_func.params = (void *)&par;

    s = gsl_multimin_fminimizer_alloc (T, n);
    gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

    size_t iterations;
    double tolerance;
    int status;
    
    iterations=0;
    tolerance=1e-5;

    const size_t itermax=3000;
    double size;

    do
    {
        iterations++;
        status = gsl_multimin_fminimizer_iterate (s);

        if (status) {
            break;
        }

        size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, tolerance);

        if (status == GSL_SUCCESS) {
            printf ("converged to minimum at\n");
            std::cout << "iterations = " << iterations << " f() = " 
            << s->fval << " size = " << size << std::endl;
        }
    } while (status == GSL_CONTINUE && iterations < itermax);

    tolerance = size;


    x = gsl_multimin_fminimizer_x(s);

    // std::cout << " iterations taken = " << iterations <<
    // " tol achieved = " << tolerance <<
    // " status = " << status << std::endl;

    // std::cout << "x vector is" << std::endl;
    // for (int i=0; i<int(n); i++) {
    //     std::cout << gsl_vector_get(x, i) << "  ";
    // }
    // std::cout << std::endl;
    // std::cout << "step size vector is" << std::endl;
    // for (int i=0; i<int(n); i++) {
    //     std::cout << gsl_vector_get(ss, i) << "  ";
    // }
    // std::cout << std::endl;

    gsl_vector_free (x);
    gsl_vector_free (ss);
    // Code crashes here.
    gsl_multimin_fminimizer_free (s);

    return(0);
}