# include <cppad/cppad.hpp>
# include <omp.h>
# include <cppad/ipopt/solve.hpp>



using namespace std;
# define NUMBER_THREADS  4

namespace {
    using CppAD::AD;

    class FG_eval {
    public:
        typedef CPPAD_TESTVECTOR(AD<double>) ADvector;
        void operator()(CPPAD_TESTVECTOR(AD<double>)& fg, const CPPAD_TESTVECTOR(AD<double>)& x)
        {
            assert(fg.size() == 3);
            assert(x.size() == 4);

            // Fortran style indexing
            AD<double> x1 = x[0];
            AD<double> x2 = x[1];
            AD<double> x3 = x[2];
            AD<double> x4 = x[3];
            // f(x)
            fg[0] = x1 * x4 * (x1 + x2 + x3) + x3;
            // g_1 (x)
            fg[1] = x1 * x2 * x3 * x4;
            // g_2 (x)
            fg[2] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
            //
            return;
        }
    };
    // structure with problem specific information
    typedef struct {
        // function argument (worker input)
        double          x;
        // This structure would also have return information in it,
        // but this example only returns the ok flag
    } problem_specific;
    // =====================================================================
    // General purpose code you can copy to your application
    // =====================================================================
    using CppAD::thread_alloc;
    // ------------------------------------------------------------------
    // used to inform CppAD when we are in parallel execution mode
    bool in_parallel(void)
    {
        return omp_in_parallel() != 0;
    }
    // ------------------------------------------------------------------
    // used to inform CppAD of the current thread number
    size_t thread_number(void)
    {
        return static_cast<size_t>(omp_get_thread_num());
    }
    // ------------------------------------------------------------------
    // structure with information for one thread
    typedef struct {
        // false if an error occurs, true otherwise (worker output)
        bool               ok;
    } thread_one_t;
    // vector with information for all threads
    thread_one_t thread_all_[NUMBER_THREADS];
    // ------------------------------------------------------------------
    // function that calls all the workers
    template <class Type>
    Type Poly(const CPPAD_TESTVECTOR(double)& a, const Type& x)
    {
        size_t k = a.size();
        Type y = 0.;  // initialize summation
        Type x_i = 1.;  // initialize x^i
        for (size_t i = 0; i < k; i++)
        {
            y += a[i] * x_i;  // y   = y + a_i * x^i
            x_i *= x;           // x_i = x_i * x
        }
        return y;
    }
    bool worker(problem_specific* info);
    bool run_all_workers(size_t num_threads, problem_specific* info_all[])
    {
        bool ok = true;

        // initialize thread_all_
        int thread_num, int_num_threads = int(num_threads);
        for (thread_num = 0; thread_num < int_num_threads; thread_num++)
        {  // initialize as false to make sure gets called for all threads
            thread_all_[thread_num].ok = false;
        }

        // turn off dynamic thread adjustment
        omp_set_dynamic(0);

        // set the number of OpenMP threads
        omp_set_num_threads(int_num_threads);

        // setup for using CppAD::AD<double> in parallel
        thread_alloc::parallel_setup(
            num_threads, in_parallel, thread_number
        );
        thread_alloc::hold_memory(true);
        CppAD::parallel_ad<double>();

        // execute worker in parallel
# pragma omp parallel for
        for (thread_num = 0; thread_num < int_num_threads; thread_num++)
            thread_all_[thread_num].ok = worker(info_all[thread_num]);
        // end omp parallel for

              // set the number of OpenMP threads to one
        omp_set_num_threads(1);

        // now inform CppAD that there is only one thread
        thread_alloc::parallel_setup(1, nullptr, nullptr);
        thread_alloc::hold_memory(false);
        CppAD::parallel_ad<double>();

        // check to ok flag returned by during calls to work by other threads
        for (thread_num = 1; thread_num < int_num_threads; thread_num++)
            ok &= thread_all_[thread_num].ok;

        return ok;
    }
    // =====================================================================
    // End of General purpose code
    // =====================================================================
    // function that does the work for one thread
    bool worker(problem_specific* info)
    {
        bool ok = true;
        size_t i;
        using CppAD::AD;
        using CppAD::vector;

        // number of independent variables (domain dimension for f and g)
        size_t nx = 4;
        // number of constraints (range dimension for g)
        size_t ng = 2;
        // initial value of the independent variables
        CPPAD_TESTVECTOR(double) xi(nx);
        xi[0] = 1.0;
        xi[1] = 5.0;
        xi[2] = 5.0;
        xi[3] = 1.0;
        // lower and upper limits for x
        CPPAD_TESTVECTOR(double) xl(nx), xu(nx);
        for (i = 0; i < nx; i++)
        {
            xl[i] = 1.0;
            xu[i] = 5.0;
        }
        // lower and upper limits for g
        CPPAD_TESTVECTOR(double) gl(ng), gu(ng);
        gl[0] = 25.0;     gu[0] = 1.0e19;
        gl[1] = 40.0;     gu[1] = 40.0;

        // object that computes objective and constraints
        FG_eval fg_eval;

        // options
        std::string options;
        // turn off any printing
        options += "Integer print_level  0\n";
        options += "String  sb           yes\n";
        // maximum number of iterations
        options += "Integer max_iter     10\n";
        // approximate accuracy in first order necessary conditions;
        // see Mathematical Programming, Volume 106, Number 1,
        // Pages 25-57, Equation (6)
        options += "Numeric tol          1e-6\n";
        // derivative testing
        options += "String  derivative_test            second-order\n";
        // maximum amount of random pertubation; e.g.,
        // when evaluation finite diff
        options += "Numeric point_perturbation_radius  0.\n";

        // place to return solution
        CppAD::ipopt::solve_result<CPPAD_TESTVECTOR(double)> solution;

        // solve the problem
        CppAD::ipopt::solve<CPPAD_TESTVECTOR(double), FG_eval>(
            options, xi, xl, xu, gl, gu, fg_eval, solution
            );
        //
        // Check some of the solution values
        //
        ok &= solution.status == CppAD::ipopt::solve_result<CPPAD_TESTVECTOR(double)>::success;
        //
        double check_x[] = { 1.000000, 4.743000, 3.82115, 1.379408 };
        double check_zl[] = { 1.087871, 0.,       0.,      0. };
        double check_zu[] = { 0.,       0.,       0.,      0. };
        double rel_tol = 1e-6;  // relative tolerance
        double abs_tol = 1e-6;  // absolute tolerance
        for (i = 0; i < nx; i++)
        {
            ok &= CppAD::NearEqual(
                check_x[i], solution.x[i], rel_tol, abs_tol
            );
            ok &= CppAD::NearEqual(
                check_zl[i], solution.zl[i], rel_tol, abs_tol
            );
            ok &= CppAD::NearEqual(
                check_zu[i], solution.zu[i], rel_tol, abs_tol
            );
        }

        return ok;
    }
}
bool simple_ad(void)
{
    bool ok = true;
    size_t num_threads = NUMBER_THREADS;

    // Check that no memory is in use or avialable at start
    // (using thread_alloc in sequential mode)
    size_t thread_num;
    for (thread_num = 0; thread_num < num_threads; thread_num++)
    {
        ok &= thread_alloc::inuse(thread_num) == 0;
        ok &= thread_alloc::available(thread_num) == 0;
    }

    // initialize info_all
    problem_specific* info, * info_all[NUMBER_THREADS];
    for (thread_num = 0; thread_num < num_threads; thread_num++)
    {  // problem specific information
        size_t min_bytes(sizeof(info)), cap_bytes;
        void* v_ptr = thread_alloc::get_memory(min_bytes, cap_bytes);
        info = static_cast<problem_specific*>(v_ptr);
        info->x = double(thread_num) + 1.;
        info_all[thread_num] = info;
    }

    ok &= run_all_workers(num_threads, info_all);

    // go down so that free memory for other threads before memory for master
    thread_num = num_threads;
    while (thread_num--)
    {  // delete problem specific information
        void* v_ptr = static_cast<void*>(info_all[thread_num]);
        thread_alloc::return_memory(v_ptr);
        // check that there is no longer any memory inuse by this thread
        ok &= thread_alloc::inuse(thread_num) == 0;
        // return all memory being held for future use by this thread
        thread_alloc::free_available(thread_num);
    }

    return ok;
}

int main() {
    bool test = simple_ad();
    cout << test << endl;
    return 0;
}