// Copyright (C) 2010 Deltares
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License version 2 as
// published by the Free Software Foundation.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

/**
 * @file
 * @brief xxx
 * @author Dirk Schwanenberg
 * @version 1.0
 * @date 2010, 2011
 */

#ifdef _MSC_VER
#define snprintf sprintf_s
#endif

#include "rtcToolsIPOPT.h"
#include "piDiagInterface.h"
#include "IpIpoptApplication.hpp"
#include "IpIpoptCalculatedQuantities.hpp"
#include "IpSolveStatistics.hpp"

const string RTCTOOLSIPOPT_CODE = "RTCTools.runtime.rtcToolsIPOPT";

#include <cstdio>
#include <stdexcept>

using namespace Ipopt;

// constructors
rtcToolsIPOPT::rtcToolsIPOPT(rtcToolsRuntime *tool, string workDir)
{
	this->tool = tool;
	this->workDir = workDir;
	this->parCustom = false;
	needsInitialization = true;
}

rtcToolsIPOPT::rtcToolsIPOPT(rtcToolsRuntime *tool, rtcRuntimeConfigSettings::IPOPT par, string workDir)
{
	this->tool = tool;
	this->workDir = workDir;
	this->parCustom = true;
	this->par = par;
	needsInitialization = true;
}

// returns the size of the problem
bool rtcToolsIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
							     Index& nnz_h_lag, IndexStyleEnum& index_style)
{
	// dimensions of optimization problem
	n = _n = tool->getN();
	m = _m = tool->getM();
	nnz_jac_g = _nnz = tool->getNNZ();

	// no updated hessian
	nnz_h_lag = 0;

	// use the C style indexing (0-based)
	index_style = TNLP::C_STYLE;

	return true;
}

// returns the variable bounds
bool rtcToolsIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u,
								    Index m, Number* g_l, Number* g_u)
{
	// here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
	// If desired, we could assert to make sure they are what we think they are.
	assert(_n == n);
	assert(_m == m);

	tool->get_bounds_info(n, (double *)x_l, (double *)x_u,
		                  m, (double *)g_l, (double *)g_u);

	return true;
}

bool rtcToolsIPOPT::get_var_con_metadata(Index n,
						  StringMetaDataMapType& var_string_md,
                          IntegerMetaDataMapType& var_integer_md,
                          NumericMetaDataMapType& var_numeric_md,
                          Index m,
                          StringMetaDataMapType& con_string_md,
                          IntegerMetaDataMapType& con_integer_md,
                          NumericMetaDataMapType& con_numeric_md)
{
	if (n>0) {
		std::vector<std::string>& varnamesvec(var_string_md["idx_names"]);
		varnamesvec.resize(n);

		vector<metadata> var(n);
		tool->getN_metadata(var);
		for (int i=0; i<n; i++) {
			varnamesvec[i] = var[i].label;
		}
	}

	if (m>0) {
		std::vector<std::string>& consnamesvec(con_string_md["idx_names"]);
		consnamesvec.resize(m);

		vector<metadata> con(m);
		tool->getM_metadata(con);
		for (int i=0; i<m; i++) {
			consnamesvec[i] = con[i].label;
		}
	}

	return true;
}

// returns the initial point for the problem
bool rtcToolsIPOPT::get_starting_point(Index n, bool init_x, Number* x,
									   bool init_z, Number* z_L, Number* z_U,
									   Index m, bool init_lambda,
									   Number* lambda)
{
	// Here, we assume we only have starting values for x, if you code
	// your own NLP, you can provide starting values for the dual variables
	// if you wish
	assert(_n == n);
	assert(init_x == true);
	assert(init_z == false);
	assert(_m == m);
	assert(init_lambda == false);

	// initialize to the given starting point
	for (Index i=0; i<n; i++) x[i] = x_0[i];

	return true;
}

double obj_value_dump;

// returns the value of the objective function
bool rtcToolsIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
{
	assert(_n == n);

	clock_t t = clock();
	if (new_x) {
		obj_value = tool->simulate(n, x);
	} else {
		obj_value = obj_value_dump;
	}
	t_eval_f += clock()-t;

	return true;
}

// return the gradient of the objective function grad_{x} f(x)
bool rtcToolsIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
{
	assert(_n == n);

	clock_t t = clock();
	if (new_x) obj_value_dump = tool->simulate(n, x);
	tool->eval_grad_f(n, grad_f);
	t_eval_grad_f += clock()-t;

	return true;
}

// return the value of the constraints: g(x)
bool rtcToolsIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
{
	assert(_n == n);
	assert(_m == m);

	clock_t t = clock();
	if (m>0) {
		if (new_x) obj_value_dump = tool->simulate(n, x);
		tool->eval_g(n, x, m, g);
	}
	t_eval_g += clock()-t;

	return true;
}

// return the structure or values of the jacobian
bool rtcToolsIPOPT::eval_jac_g(Index n, const Number* x, bool new_x,
							   Index m, Index nnz_jac_g, Index* iRow, Index *jCol,
							   Number* values)
{
	assert(_n == n);
	assert(_m == m);

	clock_t t = clock();
	if (m>0) {
		if (new_x) obj_value_dump = tool->simulate(n, x);
		tool->eval_jac_g(m, nnz_jac_g, iRow, jCol, values);
	}
	t_eval_jac_g += clock()-t;

	return true;
}

//return the structure or values of the hessian
bool rtcToolsIPOPT::eval_h(Index n, const Number* x, bool new_x,
						   Number obj_factor, Index m, const Number* lambda,
						   bool new_lambda, Index nele_hess, Index* iRow,
						   Index* jCol, Number* values)
{
	// approximated hessian
	return false;
}

bool rtcToolsIPOPT::get_scaling_parameters(Number& obj_scaling,
										   bool& use_x_scaling, Index n,
										   Number* x_scaling,
										   bool& use_g_scaling, Index m,
										   Number* g_scaling)
{
	// scaling in RTC-Tools is done internally, user supplied scaling in IPOPT is not supported
	return false;
}

bool rtcToolsIPOPT::intermediate_callback(AlgorithmMode mode,
										  Index iter, Number obj_value,
										  Number inf_pr, Number inf_du,
										  Number mu, Number d_norm,
										  Number regularization_size,
										  Number alpha_du, Number alpha_pr,
										  Index ls_trials,
										  const IpoptData* ip_data,
										  IpoptCalculatedQuantities* ip_cq)
{
	// convergence history
	hist_mode.push_back(mode);
	hist_iter.push_back(iter);
	hist_f.push_back(obj_value);
	hist_constr_viol.push_back(ip_cq->unscaled_curr_nlp_constraint_violation(NORM_MAX));
	hist_dual_inf.push_back(inf_du);
	hist_mu.push_back(mu);

	// test implementation to terminate if the constraints are fulfilled and a certain number of iterations is reached
	// termination only in regular mode, not in restoration mode
	if (iter > this->par.max_iter_user &&
		ip_cq->unscaled_curr_nlp_constraint_violation(NORM_MAX) < this->par.constr_viol_tol_user &&
		inf_du < this->par.dual_inf_tol_user &&
		mode == RegularMode) {

		char buffer[200];
		snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: user-defined criteria for constr_viol_tol (%E) and dual_inf_tol (%E) met after %d>=%d iterations",
			this->par.constr_viol_tol_user, this->par.dual_inf_tol_user, iter, this->par.max_iter_user);
		piDiagInterface::addLine(3, string(buffer));

		return false;
	}

	return true;
}

void rtcToolsIPOPT::finalize_solution(SolverReturn status,
								      Index n, const Number* x, const Number* z_L, const Number* z_U,
									  Index m, const Number* g, const Number* lambda,
									  Number obj_value,
									  const IpoptData* ip_data,
									  IpoptCalculatedQuantities* ip_cq)
{
	// ensure our timeseries are consistent with the optimal result 'x'
	tool->simulate(n, x);
}

int rtcToolsIPOPT::optimize()
{
	if (needsInitialization) {
		initialize();
		needsInitialization = false;
	}

	// gets starting point and fill NaNs with zero
	int n = tool->getN();
	x_0 = new double[n];
	tool->getInput(n, x_0);
	for (int i=0; i<n; i++) {
		if (x_0[i]!=x_0[i]) x_0[i] = 0.0;
	}

	hist_mode = vector<AlgorithmMode>();
	hist_iter = vector<Index>();
	hist_f = vector<Number>();
	hist_constr_viol = vector<Number>();
	hist_dual_inf = vector<Number>();
	hist_mu = vector<Number>();

	t_eval_f = 0;
	t_eval_grad_f = 0;
	t_eval_g = 0;
	t_eval_jac_g = 0;

	ApplicationReturnStatus status = app->OptimizeTNLP(this);

	// timing info
	char buffer[200];
	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (eval_f) = %d ms", (int)t_eval_f);
	piDiagInterface::addLine(4, string(buffer));
	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (eval_grad_f) = %d ms", (int)t_eval_grad_f);
	piDiagInterface::addLine(4, string(buffer));
	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (eval_g) = %d ms", (int)t_eval_g);
	piDiagInterface::addLine(4, string(buffer));
	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (eval_jac_g) = %d ms", (int)t_eval_jac_g);
	piDiagInterface::addLine(4, string(buffer));

	// process status message -----------------------------------------------------
	bool success = false;

	// sucessful execution
	if (status==Solve_Succeeded) {
		piDiagInterface::addLine(3, "rtcToolsIPOPT: optimization stops with Solve_Succeeded", RTCTOOLSIPOPT_CODE);
		success = true;
	} else if (status==Solved_To_Acceptable_Level) {
		piDiagInterface::addLine(3, "rtcToolsIPOPT: optimization stops with Solved_To_Acceptable_Level", RTCTOOLSIPOPT_CODE);
		success = true;
	} else if (status==User_Requested_Stop) {
		piDiagInterface::addLine(3, "rtcToolsIPOPT: optimization stops with User_Requested_Stop", RTCTOOLSIPOPT_CODE);
		success = true;
	}

	// termination, but no optimum
	else if (status==Maximum_Iterations_Exceeded) piDiagInterface::addLine(2, "rtcToolsIPOPT: optimization stops with Maximum_Iterations_Exceeded", RTCTOOLSIPOPT_CODE);
	else if (status==Maximum_CpuTime_Exceeded) piDiagInterface::addLine(2, "rtcToolsIPOPT: optimization stops with Maximum_CpuTime_Exceeded", RTCTOOLSIPOPT_CODE);
    else if (status==Search_Direction_Becomes_Too_Small) piDiagInterface::addLine(2, "rtcToolsIPOPT: optimization stops with Search_Direction_Becomes_Too_Small", RTCTOOLSIPOPT_CODE);
    else if (status==Feasible_Point_Found) piDiagInterface::addLine(2, "rtcToolsIPOPT: optimization stops with Feasible_Point_Found", RTCTOOLSIPOPT_CODE);
    else if (status==Diverging_Iterates) piDiagInterface::addLine(2, "rtcToolsIPOPT: optimization stops with Diverging_Iterates", RTCTOOLSIPOPT_CODE);

	// error during execution
	else if (status==Infeasible_Problem_Detected) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Infeasible_Problem_Detected", RTCTOOLSIPOPT_CODE);
    else if (status==Restoration_Failed) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Restoration_Failed. Unable to reduce infeasibility.", RTCTOOLSIPOPT_CODE);
    else if (status==Error_In_Step_Computation) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Error_In_Step_Computation", RTCTOOLSIPOPT_CODE);
    else if (status==Invalid_Number_Detected) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Invalid_Number_Detected", RTCTOOLSIPOPT_CODE);
    else if (status==Not_Enough_Degrees_Of_Freedom) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Not_Enough_Degrees_Of_Freedom", RTCTOOLSIPOPT_CODE);
    else if (status==Invalid_Option) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Invalid_Option", RTCTOOLSIPOPT_CODE);
    else if (status==Invalid_Problem_Definition) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Invalid_Problem_Definition", RTCTOOLSIPOPT_CODE);
    else if (status==Unrecoverable_Exception) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Unrecoverable_Exception", RTCTOOLSIPOPT_CODE);
    else if (status==NonIpopt_Exception_Thrown) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with NonIpopt_Exception_Thrown", RTCTOOLSIPOPT_CODE);
    else if (status==Insufficient_Memory) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Insufficient_Memory", RTCTOOLSIPOPT_CODE);
    else if (status==Internal_Error) piDiagInterface::addLine(1, "rtcToolsIPOPT: optimization stops with Internal_Error", RTCTOOLSIPOPT_CODE);

	else {
		stringstream ss;
		ss << "undocumented IPOPT status message " << status << " is not implemented in RTC-Tools, check the file 'ipopt.out' for details!";
		piDiagInterface::addLine(2, ss.str(), RTCTOOLSIPOPT_CODE);
	}

	if (par.treat_unsuccess_as_error && success==false) {
		static char buffer[120];
		snprintf(buffer, sizeof(buffer), "unsuccessful execution of IPOPT. status: %d", (int)status);
		throw runtime_error(buffer);
	}

	// delete array with starting point
	delete [] x_0;

	return status;
}

void rtcToolsIPOPT::initialize()
{
	app = IpoptApplicationFactory();

	// ouput
	app->Options()->SetIntegerValue("print_level", par.print_level);
	app->Options()->SetStringValue("print_user_options", par.print_user_options);
	app->Options()->SetStringValue("print_options_documentation", par.print_options_documentation);
	app->Options()->SetStringValue("print_timing_statistics", par.print_timing_statistics);
	app->Options()->SetIntegerValue("file_print_level", par.file_print_level);
	app->Options()->SetStringValue("output_file", utils::getAbsoluteFilename(workDir, "ipopt.out"));

	// termination
	app->Options()->SetNumericValue("tol", par.tol);
	app->Options()->SetIntegerValue("max_iter", par.max_iter);
	app->Options()->SetNumericValue("max_cpu_time", par.max_cpu_time);
	app->Options()->SetNumericValue("dual_inf_tol", par.dual_inf_tol);
	app->Options()->SetNumericValue("constr_viol_tol", par.constr_viol_tol);
	app->Options()->SetNumericValue("compl_inf_tol", par.compl_inf_tol);
	app->Options()->SetNumericValue("acceptable_tol", par.acceptable_tol);
	app->Options()->SetIntegerValue("acceptable_iter", par.acceptable_iter);
	app->Options()->SetNumericValue("acceptable_constr_viol_tol", par.acceptable_constr_viol_tol);
	app->Options()->SetNumericValue("acceptable_dual_inf_tol", par.acceptable_dual_inf_tol);
	app->Options()->SetNumericValue("acceptable_compl_inf_tol", par.acceptable_compl_inf_tol);
	app->Options()->SetNumericValue("acceptable_obj_change_tol", par.acceptable_obj_change_tol);
	app->Options()->SetNumericValue("diverging_iterates_tol", par.diverging_iterates_tol);

	// nlp scaling
	app->Options()->SetNumericValue("obj_scaling_factor", par.obj_scaling_factor);
	app->Options()->SetStringValue("nlp_scaling_method", par.nlp_scaling_method);
	app->Options()->SetNumericValue("nlp_scaling_max_gradient", par.nlp_scaling_max_gradient);
	app->Options()->SetNumericValue("nlp_scaling_min_value", par.nlp_scaling_min_value);

	// nlp
	app->Options()->SetNumericValue("bound_relax_factor", par.bound_relax_factor);
	app->Options()->SetStringValue("honor_original_bounds", par.honor_original_bounds);
	app->Options()->SetStringValue("check_derivatives_for_naninf", par.check_derivatives_for_naninf);
	app->Options()->SetStringValue("fixed_variable_treatment", par.fixed_variable_treatment);
	app->Options()->SetStringValue("jac_c_constant", par.jac_c_constant);
	app->Options()->SetStringValue("jac_d_constant", par.jac_d_constant);

	// derivative checker
	app->Options()->SetStringValue("derivative_test", par.derivative_test);
	app->Options()->SetNumericValue("derivative_test_perturbation", par.derivative_test_perturbation);
	app->Options()->SetNumericValue("derivative_test_tol", par.derivative_test_tol);
	app->Options()->SetStringValue("derivative_test_print_all", par.derivative_test_print_all);
	app->Options()->SetStringValue("jacobian_approximation", par.jacobian_approximation);
	app->Options()->SetNumericValue("findiff_perturbation", par.findiff_perturbation);

	app->Options()->SetStringValue("hessian_approximation", "limited-memory");
	app->Options()->SetIntegerValue("limited_memory_max_history", par.limited_memory_max_history);
	app->Options()->SetIntegerValue("limited_memory_max_skipping", par.limited_memory_max_skipping);

	// linear equation solver
	app->Options()->SetStringValue("linear_solver", par.linear_solver);
	app->Options()->SetStringValue("linear_system_scaling", par.linear_system_scaling);
	app->Options()->SetStringValue("linear_scaling_on_demand", par.linear_scaling_on_demand);
	app->Options()->SetIntegerValue("min_refinement_steps", par.min_refinement_steps);
    app->Options()->SetIntegerValue("max_refinement_steps", par.max_refinement_steps);

	// hessian perturbation
	app->Options()->SetNumericValue("min_hessian_perturbation", par.min_hessian_perturbation);
    app->Options()->SetNumericValue("max_hessian_perturbation", par.max_hessian_perturbation);
	app->Options()->SetNumericValue("first_hessian_perturbation", par.first_hessian_perturbation);
    app->Options()->SetNumericValue("perturb_inc_fact_first", par.perturb_inc_fact_first);
	app->Options()->SetNumericValue("perturb_inc_fact", par.perturb_inc_fact);
    app->Options()->SetNumericValue("perturb_dec_fact", par.perturb_dec_fact);
    app->Options()->SetNumericValue("jacobian_regularization_value", par.jacobian_regularization_value);

    // restoration phase
	app->Options()->SetStringValue("expect_infeasible_problem", par.expect_infeasible_problem);
	app->Options()->SetNumericValue("expect_infeasible_problem_ctol", par.expect_infeasible_problem_ctol);
	app->Options()->SetNumericValue("expect_infeasible_problem_ytol", par.expect_infeasible_problem_ytol);
	app->Options()->SetStringValue("start_with_resto", par.start_with_resto);
	app->Options()->SetNumericValue("soft_resto_pderror_reduction_factor", par.soft_resto_pderror_reduction_factor);
	app->Options()->SetNumericValue("required_infeasibility_reduction", par.required_infeasibility_reduction);
    app->Options()->SetNumericValue("bound_mult_reset_threshold", par.bound_mult_reset_threshold);
	app->Options()->SetNumericValue("constr_mult_reset_threshold", par.constr_mult_reset_threshold);
	app->Options()->SetStringValue("evaluate_orig_obj_at_resto_trial", par.evaluate_orig_obj_at_resto_trial);

	// multiplier updates
	app->Options()->SetStringValue("alpha_for_y", par.alpha_for_y);
	app->Options()->SetNumericValue("alpha_for_y_tol", par.alpha_for_y_tol);
	app->Options()->SetStringValue("recalc_y", par.recalc_y);
	app->Options()->SetNumericValue("recalc_y_feas_tol", par.recalc_y_feas_tol);

    // line search
	app->Options()->SetIntegerValue("max_soc", par.max_soc);
	app->Options()->SetIntegerValue("watchdog_shortened_iter_trigger", par.watchdog_shortened_iter_trigger);
	app->Options()->SetIntegerValue("watchdog_trial_iter_max", par.watchdog_trial_iter_max);
	app->Options()->SetStringValue("accept_every_trial_step", par.accept_every_trial_step);

	// initialization
	app->Options()->SetNumericValue("bound_frac", par.bound_frac);
	app->Options()->SetNumericValue("bound_push", par.bound_push);
	app->Options()->SetNumericValue("slack_bound_frac", par.slack_bound_frac);
	app->Options()->SetNumericValue("slack_bound_push", par.slack_bound_push);
	app->Options()->SetNumericValue("bound_mult_init_val", par.bound_mult_init_val);
	app->Options()->SetNumericValue("constr_mult_init_max", par.constr_mult_init_max);
    app->Options()->SetStringValue("bound_mult_init_method", par.bound_mult_init_method);

	// barrier parameter
	app->Options()->SetStringValue("mehrotra_algorithm", par.mehrotra_algorithm);
	app->Options()->SetStringValue("mu_strategy", par.mu_strategy);
	app->Options()->SetStringValue("mu_oracle", par.mu_oracle);
	app->Options()->SetIntegerValue("quality_function_max_section_steps", par.quality_function_max_section_steps);
	app->Options()->SetStringValue("fixed_mu_oracle", par.fixed_mu_oracle);
	app->Options()->SetStringValue("adaptive_mu_globalization", par.adaptive_mu_globalization);
	app->Options()->SetNumericValue("mu_init", par.mu_init);
	app->Options()->SetNumericValue("mu_max_fact", par.mu_max_fact);
	app->Options()->SetNumericValue("mu_max", par.mu_max);
	app->Options()->SetNumericValue("mu_min", par.mu_min);
	app->Options()->SetNumericValue("mu_target", par.mu_target);
	app->Options()->SetNumericValue("barrier_tol_factor", par.barrier_tol_factor);
	app->Options()->SetNumericValue("mu_linear_decrease_factor", par.mu_linear_decrease_factor);
	app->Options()->SetNumericValue("mu_superlinear_decrease_power", par.mu_superlinear_decrease_power);

	// warm start
	app->Options()->SetStringValue("warm_start_init_point", par.warm_start_init_point);
	app->Options()->SetNumericValue("warm_start_bound_frac", par.warm_start_bound_frac);
	app->Options()->SetNumericValue("warm_start_bound_push", par.warm_start_bound_push);
	app->Options()->SetNumericValue("warm_start_slack_bound_frac", par.warm_start_slack_bound_frac);
	app->Options()->SetNumericValue("warm_start_slack_bound_push", par.warm_start_slack_bound_push);
	app->Options()->SetNumericValue("warm_start_mult_bound_push", par.warm_start_mult_bound_push);
	app->Options()->SetNumericValue("warm_start_mult_init_max", par.warm_start_mult_init_max);

    // MA27 Linear Solver
	if (par.ma27options) {
		app->Options()->SetNumericValue("ma27_pivtol", par.ma27_pivtol);
		app->Options()->SetNumericValue("ma27_pivtolmax", par.ma27_pivtolmax);
		app->Options()->SetNumericValue("ma27_liw_init_factor", par.ma27_liw_init_factor);
		app->Options()->SetNumericValue("ma27_la_init_factor", par.ma27_meminc_factor);
		app->Options()->SetNumericValue("ma27_meminc_factor", par.ma27_meminc_factor);
	}

	// MA57 Linear Solver
	if (par.ma57options) {
		app->Options()->SetNumericValue("ma57_pivtol", par.ma57_pivtol);
		app->Options()->SetNumericValue("ma57_pivtolmax", par.ma57_pivtolmax);
		app->Options()->SetNumericValue("ma57_pre_alloc", par.ma57_pre_alloc);
		app->Options()->SetIntegerValue("ma57_pivot_order", par.ma57_pivot_order);
		app->Options()->SetStringValue("ma57_automatic_scaling", par.ma57_automatic_scaling);
		app->Options()->SetIntegerValue("ma57_block_size", par.ma57_block_size);
		app->Options()->SetIntegerValue("ma57_node_amalgamation", par.ma57_node_amalgamation);
		app->Options()->SetIntegerValue("ma57_small_pivot_flag", par.ma57_small_pivot_flag);
	}

	// MA77 Linear Solver
	if (par.ma77options) {
		app->Options()->SetIntegerValue("ma77_print_level", par.ma77_print_level);
		app->Options()->SetIntegerValue("ma77_buffer_lpage", par.ma77_buffer_lpage);
		app->Options()->SetIntegerValue("ma77_buffer_npage", par.ma77_buffer_npage);
		app->Options()->SetIntegerValue("ma77_file_size", par.ma77_file_size);
		app->Options()->SetIntegerValue("ma77_maxstore", par.ma77_file_size);
		app->Options()->SetIntegerValue("ma77_nemin", par.ma77_nemin);
		app->Options()->SetStringValue("ma77_order", par.ma77_order);
		app->Options()->SetNumericValue("ma77_small", par.ma77_small);
		app->Options()->SetNumericValue("ma77_static", par.ma77_static);
		app->Options()->SetNumericValue("ma77_u", par.ma77_u);
		app->Options()->SetNumericValue("ma77_umax", par.ma77_umax);
	}

	// MA86 Linear Solver
	if (par.ma86options) {
		app->Options()->SetIntegerValue("ma86_print_level", par.ma86_print_level);
		app->Options()->SetIntegerValue("ma86_nemin", par.ma86_nemin);
		app->Options()->SetStringValue("ma86_order", par.ma86_order);
		app->Options()->SetStringValue("ma86_scaling", par.ma86_scaling);
		app->Options()->SetNumericValue("ma86_small", par.ma86_small);
		app->Options()->SetNumericValue("ma86_static", par.ma86_static);
		app->Options()->SetNumericValue("ma86_u", par.ma86_u);
		app->Options()->SetNumericValue("ma86_umax", par.ma86_umax);
	}

	// MA97 Linear Solver
	if (par.ma97options) {
		app->Options()->SetIntegerValue("ma97_print_level", par.ma97_print_level);
		app->Options()->SetIntegerValue("ma97_nemin", par.ma97_nemin);
		app->Options()->SetStringValue("ma97_order", par.ma97_order);
		app->Options()->SetStringValue("ma97_scaling", par.ma97_scaling);
		app->Options()->SetStringValue("ma97_scaling1", par.ma97_scaling1);
		app->Options()->SetStringValue("ma97_scaling2", par.ma97_scaling2);
		app->Options()->SetStringValue("ma97_scaling3", par.ma97_scaling3);
		app->Options()->SetNumericValue("ma97_small", par.ma97_small);
		app->Options()->SetStringValue("ma97_solve_blas3", par.ma97_solve_blas3);
		app->Options()->SetStringValue("ma97_switch1", par.ma97_switch1);
		app->Options()->SetStringValue("ma97_switch2", par.ma97_switch2);
		app->Options()->SetStringValue("ma97_switch3", par.ma97_switch3);
		app->Options()->SetNumericValue("ma97_u", par.ma97_u);
		app->Options()->SetNumericValue("ma97_umax", par.ma97_umax);
	}

	// MUMPS Solver
	if (par.mumpsoptions) {
		app->Options()->SetNumericValue("mumps_pivtol", par.mumps_pivtol);
		app->Options()->SetNumericValue("mumps_pivtolmax", par.mumps_pivtolmax);
		app->Options()->SetIntegerValue("mumps_mem_percent", par.mumps_mem_percent);
		app->Options()->SetIntegerValue("mumps_permuting_scaling", par.mumps_permuting_scaling);
		app->Options()->SetIntegerValue("mumps_pivot_order", par.mumps_pivot_order);
		app->Options()->SetIntegerValue("mumps_scaling", par.mumps_scaling);
	}

	// Pardiso Linear Solver
	if (par.pardisooptions) {
		app->Options()->SetStringValue("pardiso_matching_strategy", par.pardiso_matching_strategy);
		app->Options()->SetIntegerValue("pardiso_max_iterative_refinement_steps" ,par.pardiso_max_iterative_refinement_steps);
		app->Options()->SetIntegerValue("pardiso_msglvl", par.pardiso_msglvl);
		app->Options()->SetStringValue("pardiso_order", par.pardiso_order);
	}

	// Intialize the IpoptApplication and process the options
	Ipopt::ApplicationReturnStatus status = app->Initialize();
	if (status != Ipopt::Solve_Succeeded) {
		piDiagInterface::addLine(1, "int main(int argc, char *argv[]) - Error during initialization of Ipopt!", RTCTOOLSIPOPT_CODE);
		throw runtime_error("int main(int argc, char *argv[]) - Error during initialization of Ipopt!");
	}
}

void rtcToolsIPOPT::write(string workDir)
{
	// csv output of convergence history
	string filename = utils::getAbsoluteFilename(workDir, "convergence_ipopt.csv");
	ofstream csvFile(filename.c_str(), ios::out | ios::trunc);

	csvFile << "iter,mode,f,constr_viol,dual_inf,mu" << endl;

	for (int i=0; i<(int)hist_f.size(); i++) {
		csvFile << hist_iter[i] << ",";
		if (hist_mode[i]==RegularMode) {
			csvFile << "RegularMode";
		} else if (hist_mode[i]==RestorationPhaseMode) {
			csvFile << "RestorationPhaseMode";
		}
		csvFile << "," << hist_f[i] << "," << hist_constr_viol[i] << "," << hist_dual_inf[i] << "," << hist_mu[i] << endl;
	}

	csvFile.close();
}