// Copyright (C) 2012 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 rtcToolsRuntime.cpp
 * @author Camiel van Breugelen, Dirk Schwanenberg
 * @version 1.0
 * @date 2010
 */

#include "rtcToolsRuntime.h"
#include "rtcToolsGAMS_OPT.h"
#include "rtcToolsIPOPT.h"
#include "rtcToolsSA.h"
#ifdef HAVE_SNOPT
#include "rtcToolsSNOPT.h"
#endif
#include "postprocess/postprocessing.h"
#include "piDiagInterface.h"
#include <utilities/utils.h>
#include "boost/filesystem.hpp"
#include "boost/date_time/posix_time/posix_time.hpp"
#include "dataBinding/pi_run.hxx"
#include "timeseries/stringContainer.h"
#include <cmath>
#include <omp.h>


#ifdef _MSC_VER
#define snprintf sprintf_s
#define strncpy strncpy_s
#endif

using namespace std;
using namespace rtctools::postprocess;
using namespace fews::PI;
using namespace boost::posix_time;

const string RTCTOOLSRUNTIME_CODE = "RTCTools.runtime.rtcToolsRuntime";


int rtcToolsRuntime::main(int argc, const char * const argv[], versionInfo info)
{
	rtcToolsRuntime* runtime;

	clock_t cpu1 = clock();

	// --- startup ----------------------------------------------------
	try {
		// runtime constructor and preparation
		boost::filesystem::path defaultDir(".");
		char schemaDir[200] = "";
		char workDir[200] = "";

		strncpy(schemaDir, defaultDir.string().c_str(), sizeof(schemaDir));
		strncpy(workDir, defaultDir.string().c_str(), sizeof(workDir));

		// try to get path from calling instance
		if (argc>0) {
		    boost::filesystem::path schemaPath(argv[0]);
		    schemaPath.remove_filename();
		    strncpy(schemaDir, schemaPath.string().c_str(), sizeof(schemaDir));
		}

		// scan input for optional schema and work directories
		for (int i=1; i<argc; i++) {
			if (!strcmp(argv[i], "-schemaDir")) {
				strncpy(schemaDir, argv[i+1], sizeof(schemaDir));
			} else if (!strcmp(argv[i], "-workDir")) {
				strncpy(workDir, argv[i+1], sizeof(workDir));
			}
		}

		runtime = new rtcToolsRuntime(schemaDir, workDir, info);

	} catch (exception &e) {
		cout << "int main(int argc, char *argv[]) - error during startup - " << e.what() << "\n";
		return 1; // failure during startup leads to direct termination
	} catch (...) {
		cout << "int main(int argc, char *argv[]) - unknown error during startup" << "\n";;
		return 1;
	}

	// set optional number of maximum threads for OpenMP parallelization
	if (runtime->runtimeSettings.nThread>0) {
		omp_set_num_threads(runtime->runtimeSettings.nThread);
	}

	runtime->cpu1 = cpu1;
	runtime->cpu2 = clock();

	// --- execution ----------------------------------------------------
	int status = 0;
	try {
		// execution
		runtime->execute();

	} catch (exception &e) {
		cout << "int main(int argc, char *argv[]) - error during execution - " << e.what() << "\n";
		status = 1; // sets flag after failure in execution and tries to write out data
	} catch (...) {
		cout << "int main(int argc, char *argv[]) - unknown error during execution" << "\n";;
		status = 1;
	}

	// --- finalization ----------------------------------------------------
	try {
		runtime->finish();
		delete runtime;

	} catch (exception &e) {
		cout << "int main(int argc, char *argv[]) - error during finalization - " << e.what() << "\n";
		status = 1;
	} catch (...) {
		cout << "int main(int argc, char *argv[]) - unknown error during finalization" << "\n";;;
		status = 1;
	}

	clock_t t3 = clock();

	return status;
}

/**
  * constructor, destructor and related classes for initialization
  */
rtcToolsRuntime::rtcToolsRuntime(const char schemaDir[], const char workDir[], versionInfo info) : info(info)
{
	string license;

	if (info.isProprietary) {
		license = "Deltares Proprietary License";

		cout << "******************************************************************************" << endl;
		cout << "       This program contains RTC-Tools, a library for real-time control.      " << endl;
		cout << "           It is released under a proprietary Deltares license.               " << endl;
		cout << "      For more information visit http://oss.deltares.nl/web/rtc-tools/        " << endl;
		cout << "******************************************************************************" << endl;
	} else {
		license = "GPL2 License";

		cout << "******************************************************************************" << endl;
		cout << "   This program contains RTC-Tools, a library for real-time control. It is    " << endl;
		cout << "  released as open source code under the GNU General Public License 2 (GPL2). " << endl;
		cout << "      For more information visit http://oss.deltares.nl/web/rtc-tools/        " << endl;
		cout << "******************************************************************************" << endl;
	}

	eval_g_new = true;

    iEnd = 0;

	runtimeSettings.schemaDir = utils::getAbsoluteFilename(string(schemaDir));
	runtimeSettings.workDir = utils::getAbsoluteFilename(string(workDir));

	piDiagInterface::initializeRtcDiagWriter(runtimeSettings.workDir);

	// model identifier
	piDiagInterface::addLine(3, "rtcToolsRuntime: RTC-Tools, version "
		+ info.version + "." + info.revision + " (" + license + ")");
	piDiagInterface::addLine(3, "rtcToolsRuntime: model execution started at "
		+ to_simple_string(second_clock::universal_time()) + " UTC time");

	try {
		// check availability of schemas
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_diag.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_modelparameters.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_sharedtypes.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_state.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_timeseries.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcDataConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcObjectiveConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcRuntimeConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcScenarioTreeConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcSharedTypes.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcStateConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcToolsConfig.xsd"));
		assertFile(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "treeVector.xsd"));

		// read rtcRuntimeConfig.xml and other input files
		parseRuntimeConfigFile();
		parseInputFiles();

	} catch (exception &e) {
		piDiagInterface::addLine(1, "rtcToolsRuntime::rtcToolsRuntime(...) - error - " + string(e.what()));
		piDiagInterface::write();
		throw;
	} catch (...) {
		piDiagInterface::addLine(1, "rtcToolsRuntime::rtcToolsRuntime(...) - unknown error.");
		piDiagInterface::write();
		throw;
	}
}


rtcToolsRuntime::~rtcToolsRuntime(void)
{
	delete obj;
	delete schema;
	delete parInt;
	delete tsModel;
	delete[] rtcSim;
}


rtcRuntimeConfigSettings::GAMS getGAMSSettings(fews::OptimizerComplexType::GAMSType gams, parameterInterface *parInt)
{
	rtcRuntimeConfigSettings::GAMS optimizerGAMS;

	optimizerGAMS.algorithm = gams.getAlgorithmDefaultValue();
	if (gams.getAlgorithm().present()) optimizerGAMS.algorithm = gams.getAlgorithm().get();

	return optimizerGAMS;
}


rtcRuntimeConfigSettings::IPOPT getIPOPTSettings(fews::OptimizerComplexType::IPOPTType ipopt, parameterInterface *parInt, versionInfo info)
{
	rtcRuntimeConfigSettings::IPOPT optimizerIPOPT;

	// output
	IPOPTComplexType::OutputType output = ipopt.getOutput();

	optimizerIPOPT.print_level = (int)output.getPrint_levelDefaultValue();
	if (output.getPrint_level().present()) optimizerIPOPT.print_level = (int)output.getPrint_level().get();

	optimizerIPOPT.print_user_options = output.getPrint_user_optionsDefaultValue();
	if (output.getPrint_user_options().present()) optimizerIPOPT.print_user_options = output.getPrint_user_options().get();

	optimizerIPOPT.print_options_documentation = output.getPrint_options_documentationDefaultValue();
	if (output.getPrint_options_documentation().present()) optimizerIPOPT.print_options_documentation = output.getPrint_options_documentation().get();

	optimizerIPOPT.print_timing_statistics = output.getPrint_timing_statisticsDefaultValue();
	if (output.getPrint_timing_statistics().present()) optimizerIPOPT.print_timing_statistics = output.getPrint_timing_statistics().get();

	optimizerIPOPT.file_print_level = (int)output.getFile_print_levelDefaultValue();
	if (output.getFile_print_level().present()) optimizerIPOPT.file_print_level = (int)output.getFile_print_level().get();

	// termination
	IPOPTComplexType::TerminationType term = ipopt.getTermination();

	optimizerIPOPT.tol = parInt->getDblParameter(term.getTolDefaultValue());
	if (term.getTol().present()) optimizerIPOPT.tol = parInt->getDblParameter(term.getTol().get());

	optimizerIPOPT.max_iter = parInt->getIntParameter(term.getMax_iterDefaultValue());
	if (term.getMax_iter().present()) optimizerIPOPT.max_iter = parInt->getIntParameter(term.getMax_iter().get());

	optimizerIPOPT.max_cpu_time = parInt->getDblParameter(term.getMax_cpu_timeDefaultValue());
	if (term.getMax_cpu_time().present()) optimizerIPOPT.max_cpu_time = parInt->getDblParameter(term.getMax_cpu_time().get());

	optimizerIPOPT.dual_inf_tol = parInt->getDblParameter(term.getDual_inf_tolDefaultValue());
	if (term.getDual_inf_tol().present()) optimizerIPOPT.dual_inf_tol = parInt->getDblParameter(term.getDual_inf_tol().get());

	optimizerIPOPT.constr_viol_tol = parInt->getDblParameter(term.getConstr_viol_tolDefaultValue());
	if (term.getConstr_viol_tol().present()) optimizerIPOPT.constr_viol_tol = parInt->getDblParameter(term.getConstr_viol_tol().get());

	optimizerIPOPT.compl_inf_tol = parInt->getDblParameter(term.getCompl_inf_tolDefaultValue());
	if (term.getCompl_inf_tol().present()) optimizerIPOPT.compl_inf_tol = parInt->getDblParameter(term.getCompl_inf_tol().get());

	optimizerIPOPT.acceptable_tol = parInt->getDblParameter(term.getAcceptable_tolDefaultValue());
	if (term.getAcceptable_tol().present()) optimizerIPOPT.acceptable_tol = parInt->getDblParameter(term.getAcceptable_tol().get());

	optimizerIPOPT.acceptable_iter = parInt->getIntParameter(term.getAcceptable_iterDefaultValue());
	if (term.getAcceptable_iter().present()) optimizerIPOPT.acceptable_iter = parInt->getIntParameter(term.getAcceptable_iter().get());

	optimizerIPOPT.acceptable_constr_viol_tol = parInt->getDblParameter(term.getAcceptable_constr_viol_tolDefaultValue());
	if (term.getAcceptable_constr_viol_tol().present()) optimizerIPOPT.acceptable_constr_viol_tol = parInt->getDblParameter(term.getAcceptable_constr_viol_tol().get());

	optimizerIPOPT.acceptable_dual_inf_tol = parInt->getDblParameter(term.getAcceptable_dual_inf_tolDefaultValue());
	if (term.getAcceptable_dual_inf_tol().present()) optimizerIPOPT.acceptable_dual_inf_tol = parInt->getDblParameter(term.getAcceptable_dual_inf_tol().get());

	optimizerIPOPT.acceptable_compl_inf_tol = parInt->getDblParameter(term.getAcceptable_compl_inf_tolDefaultValue());
	if (term.getAcceptable_compl_inf_tol().present()) optimizerIPOPT.acceptable_compl_inf_tol = parInt->getDblParameter(term.getAcceptable_compl_inf_tol().get());

	optimizerIPOPT.acceptable_obj_change_tol = parInt->getDblParameter(term.getAcceptable_obj_change_tolDefaultValue());
	if (term.getAcceptable_obj_change_tol().present()) optimizerIPOPT.acceptable_obj_change_tol = parInt->getDblParameter(term.getAcceptable_obj_change_tol().get());

	optimizerIPOPT.diverging_iterates_tol = parInt->getDblParameter(term.getDiverging_iterates_tolDefaultValue());
	if (term.getDiverging_iterates_tol().present()) optimizerIPOPT.diverging_iterates_tol = parInt->getDblParameter(term.getDiverging_iterates_tol().get());

	optimizerIPOPT.treat_unsuccess_as_error = false;
	if (term.getTreat_unsuccess_as_error().present()) optimizerIPOPT.treat_unsuccess_as_error = parInt->getBoolParameter(term.getTreat_unsuccess_as_error().get());

	// termination (user requested)
	optimizerIPOPT.max_iter_user = 1000000;
	if (ipopt.getTerminationUser().present()) {
		optimizerIPOPT.max_iter_user = parInt->getIntParameter(ipopt.getTerminationUser().get().getMax_iter());
		optimizerIPOPT.dual_inf_tol_user = parInt->getDblParameter(ipopt.getTerminationUser().get().getDual_inf_tol());
		optimizerIPOPT.constr_viol_tol_user = parInt->getDblParameter(ipopt.getTerminationUser().get().getConstr_viol_tol());
	}

	// nlp scaling
	IPOPTComplexType::NlpScalingType scale = ipopt.getNlpScaling();

	optimizerIPOPT.obj_scaling_factor = parInt->getDblParameter(scale.getObj_scaling_factorDefaultValue());
	if (scale.getObj_scaling_factor().present()) optimizerIPOPT.obj_scaling_factor = parInt->getDblParameter(scale.getObj_scaling_factor().get());

	optimizerIPOPT.nlp_scaling_method = scale.getNlp_scaling_methodDefaultValue();
	if (scale.getNlp_scaling_method().present()) optimizerIPOPT.nlp_scaling_method = scale.getNlp_scaling_method().get();

	optimizerIPOPT.nlp_scaling_max_gradient = scale.getNlp_scaling_max_gradientDefaultValue();
	if (scale.getNlp_scaling_max_gradient().present()) optimizerIPOPT.nlp_scaling_max_gradient = scale.getNlp_scaling_max_gradient().get();

	optimizerIPOPT.nlp_scaling_min_value = scale.getNlp_scaling_min_valueDefaultValue();
	if (scale.getNlp_scaling_min_value().present()) optimizerIPOPT.nlp_scaling_min_value = scale.getNlp_scaling_min_value().get();

	// nlp
	IPOPTComplexType::NlpType nlp = ipopt.getNlp();

	optimizerIPOPT.bound_relax_factor = nlp.getBound_relax_factorDefaultValue();
	if (nlp.getBound_relax_factor().present()) optimizerIPOPT.bound_relax_factor = nlp.getBound_relax_factor().get();

	optimizerIPOPT.honor_original_bounds = nlp.getHonor_original_boundsDefaultValue();
	if (nlp.getHonor_original_bounds().present()) optimizerIPOPT.honor_original_bounds = nlp.getHonor_original_bounds().get();

	optimizerIPOPT.check_derivatives_for_naninf = nlp.getCheck_derivatives_for_naninfDefaultValue();
	if (nlp.getCheck_derivatives_for_naninf().present()) optimizerIPOPT.check_derivatives_for_naninf = nlp.getCheck_derivatives_for_naninf().get();

	optimizerIPOPT.fixed_variable_treatment = nlp.getFixed_variable_treatmentDefaultValue();
	if (nlp.getFixed_variable_treatment().present()) optimizerIPOPT.fixed_variable_treatment = nlp.getFixed_variable_treatment().get();

	optimizerIPOPT.jac_c_constant = nlp.getJac_c_constantDefaultValue();
	if (nlp.getJac_c_constant().present()) optimizerIPOPT.jac_c_constant = nlp.getJac_c_constant().get();

	optimizerIPOPT.jac_d_constant = nlp.getJac_d_constantDefaultValue();
	if (nlp.getJac_d_constant().present()) optimizerIPOPT.jac_d_constant = nlp.getJac_d_constant().get();

	// linear solver
	IPOPTComplexType::LinearSolverType ls = ipopt.getLinearSolver();

	optimizerIPOPT.linear_solver = ls.getLinear_solverDefaultValue();
	if (ls.getLinear_solver().present()) optimizerIPOPT.linear_solver = ls.getLinear_solver().get();

	optimizerIPOPT.linear_system_scaling = ls.getLinear_system_scalingDefaultValue();
	if (ls.getLinear_system_scaling().present()) optimizerIPOPT.linear_system_scaling = ls.getLinear_system_scaling().get();

	optimizerIPOPT.linear_scaling_on_demand = ls.getLinear_scaling_on_demandDefaultValue();
	if (ls.getLinear_scaling_on_demand().present()) optimizerIPOPT.linear_scaling_on_demand = ls.getLinear_scaling_on_demand().get();

	optimizerIPOPT.min_refinement_steps = ls.getMin_refinement_stepsDefaultValue();
	if (ls.getMin_refinement_steps().present()) optimizerIPOPT.min_refinement_steps = ls.getMin_refinement_steps().get();

	optimizerIPOPT.max_refinement_steps = ls.getMax_refinement_stepsDefaultValue();
	if (ls.getMax_refinement_steps().present()) optimizerIPOPT.max_refinement_steps = ls.getMax_refinement_steps().get();

	// MA27 Linear Solver
	optimizerIPOPT.ma27options = false;

	if (ls.getMa27().present()) {

		optimizerIPOPT.ma27options = true;

		optimizerIPOPT.ma27_pivtol = ls.getMa27().get().getMa27_pivtolDefaultValue();
		if (ls.getMa27().get().getMa27_pivtol().present()) optimizerIPOPT.ma27_pivtol = ls.getMa27().get().getMa27_pivtol().get();

		optimizerIPOPT.ma27_pivtolmax = ls.getMa27().get().getMa27_pivtolmaxDefaultValue();
		if (ls.getMa27().get().getMa27_pivtolmax().present()) optimizerIPOPT.ma27_pivtolmax = ls.getMa27().get().getMa27_pivtolmax().get();

		optimizerIPOPT.ma27_liw_init_factor = ls.getMa27().get().getMa27_liw_init_factorDefaultValue();
		if (ls.getMa27().get().getMa27_liw_init_factor().present()) optimizerIPOPT.ma27_liw_init_factor = ls.getMa27().get().getMa27_liw_init_factor().get();

		optimizerIPOPT.ma27_la_init_factor = ls.getMa27().get().getMa27_la_init_factorDefaultValue();
		if (ls.getMa27().get().getMa27_la_init_factor().present()) optimizerIPOPT.ma27_la_init_factor = ls.getMa27().get().getMa27_la_init_factor().get();

		optimizerIPOPT.ma27_meminc_factor = ls.getMa27().get().getMa27_meminc_factorDefaultValue();
		if (ls.getMa27().get().getMa27_meminc_factor().present()) optimizerIPOPT.ma27_meminc_factor = ls.getMa27().get().getMa27_meminc_factor().get();
	}

	// MA57 Linear Solver
	optimizerIPOPT.ma57options = false;

	if (ls.getMa57().present()) {

		optimizerIPOPT.ma57options = true;

		optimizerIPOPT.ma57_pivtol = ls.getMa57().get().getMa57_pivtolDefaultValue();
		if (ls.getMa57().get().getMa57_pivtol().present()) optimizerIPOPT.ma57_pivtol = ls.getMa57().get().getMa57_pivtol().get();

		optimizerIPOPT.ma57_pivtolmax = ls.getMa57().get().getMa57_pivtolmaxDefaultValue();
		if (ls.getMa57().get().getMa57_pivtolmax().present()) optimizerIPOPT.ma57_pivtolmax = ls.getMa57().get().getMa57_pivtolmax().get();

		optimizerIPOPT.ma57_pre_alloc = ls.getMa57().get().getMa57_pre_allocDefaultValue();
		if (ls.getMa57().get().getMa57_pre_alloc().present()) optimizerIPOPT.ma57_pre_alloc = ls.getMa57().get().getMa57_pre_alloc().get();

		optimizerIPOPT.ma57_pivot_order = ls.getMa57().get().getMa57_pivot_orderDefaultValue();
		if (ls.getMa57().get().getMa57_pivot_order().present()) optimizerIPOPT.ma57_pivot_order = ls.getMa57().get().getMa57_pivot_order().get();

		optimizerIPOPT.ma57_automatic_scaling = ls.getMa57().get().getMa57_automatic_scalingDefaultValue();
		if (ls.getMa57().get().getMa57_automatic_scaling().present()) optimizerIPOPT.ma57_automatic_scaling = ls.getMa57().get().getMa57_automatic_scaling().get();

		optimizerIPOPT.ma57_block_size = ls.getMa57().get().getMa57_block_sizeDefaultValue();
		if (ls.getMa57().get().getMa57_block_size().present()) optimizerIPOPT.ma57_block_size = ls.getMa57().get().getMa57_block_size().get();

		optimizerIPOPT.ma57_node_amalgamation = ls.getMa57().get().getMa57_node_amalgamationDefaultValue();
		if (ls.getMa57().get().getMa57_node_amalgamation().present()) optimizerIPOPT.ma57_node_amalgamation = ls.getMa57().get().getMa57_node_amalgamation().get();

		optimizerIPOPT.ma57_small_pivot_flag = ls.getMa57().get().getMa57_small_pivot_flagDefaultValue();
		if (ls.getMa57().get().getMa57_small_pivot_flag().present()) optimizerIPOPT.ma57_small_pivot_flag = ls.getMa57().get().getMa57_small_pivot_flag().get();

	}

	// MA77 Linear Solver
	optimizerIPOPT.ma77options = false;

	if (ls.getMa77().present()) {
		optimizerIPOPT.ma77options = true;

		optimizerIPOPT.ma77_print_level = ls.getMa77().get().getMa77_print_levelDefaultValue();
		if (ls.getMa77().get().getMa77_print_level().present()) optimizerIPOPT.ma77_print_level = ls.getMa77().get().getMa77_print_level().get();

	     optimizerIPOPT.ma77_buffer_lpage = ls.getMa77().get().getMa77_buffer_lpageDefaultValue();
		if (ls.getMa77().get().getMa77_buffer_lpage().present()) optimizerIPOPT.ma77_buffer_lpage = ls.getMa77().get().getMa77_buffer_lpage().get();

		 optimizerIPOPT.ma77_buffer_npage = ls.getMa77().get().getMa77_buffer_npageDefaultValue();
		if (ls.getMa77().get().getMa77_buffer_npage().present()) optimizerIPOPT.ma77_buffer_npage = ls.getMa77().get().getMa77_buffer_npage().get();

		optimizerIPOPT.ma77_file_size = ls.getMa77().get().getMa77_file_sizeDefaultValue();
		if (ls.getMa77().get().getMa77_file_size().present()) optimizerIPOPT.ma77_file_size = ls.getMa77().get().getMa77_file_size().get();

		optimizerIPOPT.ma77_maxstore = ls.getMa77().get().getMa77_maxstoreDefaultValue();
		if (ls.getMa77().get().getMa77_maxstore().present()) optimizerIPOPT.ma77_maxstore = ls.getMa77().get().getMa77_maxstore().get();

		optimizerIPOPT.ma77_nemin = ls.getMa77().get().getMa77_neminDefaultValue();
		if (ls.getMa77().get().getMa77_nemin().present()) optimizerIPOPT.ma77_nemin = ls.getMa77().get().getMa77_nemin().get();

		optimizerIPOPT.ma77_order = ls.getMa77().get().getMa77_orderDefaultValue();
		if (ls.getMa77().get().getMa77_order().present()) optimizerIPOPT.ma77_order = ls.getMa77().get().getMa77_order().get();

		optimizerIPOPT.ma77_small = ls.getMa77().get().getMa77_smallDefaultValue();
		if (ls.getMa77().get().getMa77_small().present()) optimizerIPOPT.ma77_small = ls.getMa77().get().getMa77_small().get();

		optimizerIPOPT.ma77_static = ls.getMa77().get().getMa77_staticDefaultValue();
		if (ls.getMa77().get().getMa77_static().present()) optimizerIPOPT.ma77_static = ls.getMa77().get().getMa77_static().get();

		optimizerIPOPT.ma77_u = ls.getMa77().get().getMa77_uDefaultValue();
		if (ls.getMa77().get().getMa77_u().present()) optimizerIPOPT.ma77_u = ls.getMa77().get().getMa77_u().get();

		optimizerIPOPT.ma77_umax = ls.getMa77().get().getMa77_umaxDefaultValue();
		if (ls.getMa77().get().getMa77_umax().present()) optimizerIPOPT.ma77_umax = ls.getMa77().get().getMa77_umax().get();
	}

	// MA86 Linear Solver
	optimizerIPOPT.ma86options = false;

	if (ls.getMa86().present()) {
		optimizerIPOPT.ma86options = true;

		optimizerIPOPT.ma86_print_level = ls.getMa86().get().getMa86_print_levelDefaultValue();
		if (ls.getMa86().get().getMa86_print_level().present()) optimizerIPOPT.ma86_print_level = ls.getMa86().get().getMa86_print_level().get();

		optimizerIPOPT.ma86_nemin = ls.getMa86().get().getMa86_neminDefaultValue();
		if (ls.getMa86().get().getMa86_nemin().present()) optimizerIPOPT.ma86_nemin = ls.getMa86().get().getMa86_nemin().get();

		optimizerIPOPT.ma86_order = ls.getMa86().get().getMa86_orderDefaultValue();
		if (ls.getMa86().get().getMa86_order().present()) optimizerIPOPT.ma86_order = ls.getMa86().get().getMa86_order().get();

		optimizerIPOPT.ma86_scaling = ls.getMa86().get().getMa86_scalingDefaultValue();
		if (ls.getMa86().get().getMa86_scaling().present()) optimizerIPOPT.ma86_scaling = ls.getMa86().get().getMa86_scaling().get();

		optimizerIPOPT.ma86_small = ls.getMa86().get().getMa86_smallDefaultValue();
		if (ls.getMa86().get().getMa86_small().present()) optimizerIPOPT.ma86_small = ls.getMa86().get().getMa86_small().get();

		optimizerIPOPT.ma86_static = ls.getMa86().get().getMa86_staticDefaultValue();
		if (ls.getMa86().get().getMa86_static().present()) optimizerIPOPT.ma86_static = ls.getMa86().get().getMa86_static().get();

		optimizerIPOPT.ma86_u = ls.getMa86().get().getMa86_uDefaultValue();
		if (ls.getMa86().get().getMa86_u().present()) optimizerIPOPT.ma86_u = ls.getMa86().get().getMa86_u().get();

		optimizerIPOPT.ma86_umax = ls.getMa86().get().getMa86_umaxDefaultValue();
		if (ls.getMa86().get().getMa86_umax().present()) optimizerIPOPT.ma86_umax = ls.getMa86().get().getMa86_umax().get();
	}

    // MA97 Linear Solver
	optimizerIPOPT.ma97options = false;

	if (ls.getMa97().present()) {

		optimizerIPOPT.ma97options = true;

		optimizerIPOPT.ma97_print_level = ls.getMa97().get().getMa97_print_levelDefaultValue();
		if (ls.getMa97().get().getMa97_print_level().present()) optimizerIPOPT.ma97_print_level = ls.getMa97().get().getMa97_print_level().get();

		optimizerIPOPT.ma97_nemin = ls.getMa97().get().getMa97_neminDefaultValue();
		if (ls.getMa97().get().getMa97_nemin().present()) optimizerIPOPT.ma97_nemin = ls.getMa97().get().getMa97_nemin().get();

		optimizerIPOPT.ma97_order = ls.getMa97().get().getMa97_orderDefaultValue();
		if (ls.getMa97().get().getMa97_order().present()) optimizerIPOPT.ma97_order = ls.getMa97().get().getMa97_order().get();

		optimizerIPOPT.ma97_scaling = ls.getMa97().get().getMa97_scalingDefaultValue();
		if (ls.getMa97().get().getMa97_scaling().present()) optimizerIPOPT.ma97_scaling = ls.getMa97().get().getMa97_scaling().get();

		optimizerIPOPT.ma97_scaling1 = ls.getMa97().get().getMa97_scaling1DefaultValue();
		if (ls.getMa97().get().getMa97_scaling1().present()) optimizerIPOPT.ma97_scaling1 = ls.getMa97().get().getMa97_scaling1().get();

		optimizerIPOPT.ma97_scaling2 = ls.getMa97().get().getMa97_scaling2DefaultValue();
		if (ls.getMa97().get().getMa97_scaling2().present()) optimizerIPOPT.ma97_scaling2 = ls.getMa97().get().getMa97_scaling2().get();

		optimizerIPOPT.ma97_scaling3 = ls.getMa97().get().getMa97_scaling3DefaultValue();
		if (ls.getMa97().get().getMa97_scaling3().present()) optimizerIPOPT.ma97_scaling3 = ls.getMa97().get().getMa97_scaling3().get();

		optimizerIPOPT.ma97_small = ls.getMa97().get().getMa97_smallDefaultValue();
		if (ls.getMa97().get().getMa97_small().present()) optimizerIPOPT.ma97_small = ls.getMa97().get().getMa97_small().get();

		optimizerIPOPT.ma97_solve_blas3 = ls.getMa97().get().getMa97_solve_blas3DefaultValue();
		if (ls.getMa97().get().getMa97_solve_blas3().present()) optimizerIPOPT.ma97_solve_blas3 = ls.getMa97().get().getMa97_solve_blas3().get();

		optimizerIPOPT.ma97_switch1 = ls.getMa97().get().getMa97_switch1DefaultValue();
		if (ls.getMa97().get().getMa97_switch1().present()) optimizerIPOPT.ma97_switch1 = ls.getMa97().get().getMa97_switch1().get();

		optimizerIPOPT.ma97_switch2 = ls.getMa97().get().getMa97_switch2DefaultValue();
		if (ls.getMa97().get().getMa97_switch2().present()) optimizerIPOPT.ma97_switch2 = ls.getMa97().get().getMa97_switch2().get();

		optimizerIPOPT.ma97_switch3 = ls.getMa97().get().getMa97_switch3DefaultValue();
		if (ls.getMa97().get().getMa97_switch3().present()) optimizerIPOPT.ma97_switch3 = ls.getMa97().get().getMa97_switch3().get();

		optimizerIPOPT.ma97_u = ls.getMa97().get().getMa97_uDefaultValue();
		if (ls.getMa97().get().getMa97_u().present()) optimizerIPOPT.ma97_u = ls.getMa97().get().getMa97_u().get();

		optimizerIPOPT.ma97_umax = ls.getMa97().get().getMa97_umaxDefaultValue();
		if (ls.getMa97().get().getMa97_umax().present()) optimizerIPOPT.ma97_umax = ls.getMa97().get().getMa97_umax().get();
	}

	// Mumps Linear Solver
	optimizerIPOPT.mumpsoptions = false;

	if (ls.getMumps().present()) {
		optimizerIPOPT.mumpsoptions = true;

		optimizerIPOPT.mumps_pivtol = ls.getMumps().get().getMumps_pivtolDefaultValue();
		if (ls.getMumps().get().getMumps_pivtol().present()) optimizerIPOPT.mumps_pivtol = ls.getMumps().get().getMumps_pivtol().get();

		optimizerIPOPT.mumps_pivtolmax = ls.getMumps().get().getMumps_pivtolmaxDefaultValue();
		if (ls.getMumps().get().getMumps_pivtolmax().present()) optimizerIPOPT.mumps_pivtolmax = ls.getMumps().get().getMumps_pivtolmax().get();

		optimizerIPOPT.mumps_mem_percent = ls.getMumps().get().getMumps_mem_percentDefaultValue();
		if (ls.getMumps().get().getMumps_mem_percent().present()) optimizerIPOPT.mumps_mem_percent = ls.getMumps().get().getMumps_mem_percent().get();

		optimizerIPOPT.mumps_permuting_scaling = ls.getMumps().get().getMumps_permuting_scalingDefaultValue();
		if (ls.getMumps().get().getMumps_permuting_scaling().present()) optimizerIPOPT.mumps_permuting_scaling = ls.getMumps().get().getMumps_permuting_scaling().get();

		optimizerIPOPT.mumps_pivot_order = ls.getMumps().get().getMumps_pivot_orderDefaultValue();
		if (ls.getMumps().get().getMumps_pivot_order().present()) optimizerIPOPT.mumps_pivot_order = ls.getMumps().get().getMumps_pivot_order().get();

		optimizerIPOPT.mumps_scaling = ls.getMumps().get().getMumps_scalingDefaultValue();
		if (ls.getMumps().get().getMumps_scaling().present()) optimizerIPOPT.mumps_scaling = ls.getMumps().get().getMumps_scaling().get();
	}

	// Pardiso Linear Solver
	optimizerIPOPT.pardisooptions = false;

	if (ls.getPardiso().present()) {

		optimizerIPOPT.pardisooptions = true;

		optimizerIPOPT.pardiso_matching_strategy = ls.getPardiso().get().getPardiso_matching_strategyDefaultValue();
		if (ls.getPardiso().get().getPardiso_matching_strategy().present()) optimizerIPOPT.pardiso_matching_strategy = ls.getPardiso().get().getPardiso_matching_strategy().get();

		optimizerIPOPT.pardiso_max_iterative_refinement_steps = ls.getPardiso().get().getPardiso_max_iterative_refinement_stepsDefaultValue();
		if (ls.getPardiso().get().getPardiso_max_iterative_refinement_steps().present()) optimizerIPOPT.pardiso_max_iterative_refinement_steps = ls.getPardiso().get().getPardiso_max_iterative_refinement_steps().get();

		optimizerIPOPT.pardiso_msglvl = ls.getPardiso().get().getPardiso_msglvlDefaultValue();
		if (ls.getPardiso().get().getPardiso_msglvl().present()) optimizerIPOPT.pardiso_msglvl = ls.getPardiso().get().getPardiso_msglvl().get();

		optimizerIPOPT.pardiso_order = ls.getPardiso().get().getPardiso_orderDefaultValue();
		if (ls.getPardiso().get().getPardiso_order().present()) optimizerIPOPT.pardiso_order = ls.getPardiso().get().getPardiso_order().get();
	}

	//hessian perturbation
	IPOPTComplexType::HessianPermutationType hp = ipopt.getHessianPermutation();

	optimizerIPOPT.min_hessian_perturbation = hp.getMin_hessian_perturbationDefaultValue();
	if (hp.getMin_hessian_perturbation().present()) optimizerIPOPT.min_hessian_perturbation = hp.getMin_hessian_perturbation().get();

	optimizerIPOPT.max_hessian_perturbation = hp.getMax_hessian_perturbationDefaultValue();
	if (hp.getMax_hessian_perturbation().present()) optimizerIPOPT.max_hessian_perturbation = hp.getMax_hessian_perturbation().get();

	optimizerIPOPT.first_hessian_perturbation = hp.getFirst_hessian_perturbationDefaultValue();
	if (hp.getFirst_hessian_perturbation().present()) optimizerIPOPT.first_hessian_perturbation = hp.getFirst_hessian_perturbation().get();

	optimizerIPOPT.perturb_inc_fact_first = hp.getPerturb_inc_fact_firstDefaultValue();
	if (hp.getPerturb_inc_fact_first().present()) optimizerIPOPT.perturb_inc_fact_first = hp.getPerturb_inc_fact_first().get();

    optimizerIPOPT.perturb_inc_fact = hp.getPerturb_inc_factDefaultValue();
	if (hp.getPerturb_inc_fact().present()) optimizerIPOPT.perturb_inc_fact = hp.getPerturb_inc_fact().get();

    optimizerIPOPT.perturb_dec_fact = hp.getPerturb_dec_factDefaultValue();
	if (hp.getPerturb_dec_fact().present()) optimizerIPOPT.perturb_dec_fact = hp.getPerturb_dec_fact().get();

    optimizerIPOPT.jacobian_regularization_value = hp.getJacobian_regularization_valueDefaultValue();
	if (hp.getJacobian_regularization_value().present()) optimizerIPOPT.jacobian_regularization_value = hp.getJacobian_regularization_value().get();

	//restoration phase
	IPOPTComplexType:: RestorationPhaseType rp = ipopt.getRestorationPhase();

	optimizerIPOPT.expect_infeasible_problem = rp.getExpect_infeasible_problemDefaultValue();
	if (rp.getExpect_infeasible_problem().present()) optimizerIPOPT.expect_infeasible_problem = rp.getExpect_infeasible_problem().get();

	optimizerIPOPT.expect_infeasible_problem_ctol = rp.getExpect_infeasible_problem_ctolDefaultValue();
	if (rp.getExpect_infeasible_problem_ctol().present()) optimizerIPOPT.expect_infeasible_problem_ctol = rp.getExpect_infeasible_problem_ctol().get();

	optimizerIPOPT.expect_infeasible_problem_ytol = rp.getExpect_infeasible_problem_ytolDefaultValue();
	if (rp.getExpect_infeasible_problem_ytol().present()) optimizerIPOPT.expect_infeasible_problem_ytol = rp.getExpect_infeasible_problem_ytol().get();

	optimizerIPOPT.start_with_resto = rp.getStart_with_restoDefaultValue();
	if (rp.getStart_with_resto().present()) optimizerIPOPT.start_with_resto = rp.getStart_with_resto().get();

	optimizerIPOPT.soft_resto_pderror_reduction_factor = rp.getSoft_resto_pderror_reduction_factorDefaultValue();
	if (rp.getSoft_resto_pderror_reduction_factor().present()) optimizerIPOPT.soft_resto_pderror_reduction_factor = rp.getSoft_resto_pderror_reduction_factor().get();

	optimizerIPOPT.required_infeasibility_reduction = rp.getRequired_infeasibility_reductionDefaultValue();
	if (rp.getRequired_infeasibility_reduction().present()) optimizerIPOPT.required_infeasibility_reduction = rp.getRequired_infeasibility_reduction().get();

	optimizerIPOPT.bound_mult_reset_threshold = rp.getBound_mult_reset_thresholdDefaultValue();
	if (rp.getBound_mult_reset_threshold().present()) optimizerIPOPT.bound_mult_reset_threshold = rp.getBound_mult_reset_threshold().get();

	optimizerIPOPT.constr_mult_reset_threshold = rp.getConstr_mult_reset_thresholdDefaultValue();
	if (rp.getConstr_mult_reset_threshold().present()) optimizerIPOPT.constr_mult_reset_threshold = rp.getConstr_mult_reset_threshold().get();

	optimizerIPOPT.evaluate_orig_obj_at_resto_trial = rp.getEvaluate_orig_obj_at_resto_trialDefaultValue();
	if (rp.getEvaluate_orig_obj_at_resto_trial().present()) optimizerIPOPT.evaluate_orig_obj_at_resto_trial = rp.getEvaluate_orig_obj_at_resto_trial().get();

	//multiplier updates
	IPOPTComplexType:: MultiplierUpdatesType mu = ipopt.getMultiplierUpdates();

	optimizerIPOPT.alpha_for_y = mu.getAlpha_for_yDefaultValue();
	if (mu.getAlpha_for_y().present()) optimizerIPOPT.alpha_for_y = mu.getAlpha_for_y().get();

	optimizerIPOPT.alpha_for_y_tol = mu.getAlpha_for_y_tolDefaultValue();
	if (mu.getAlpha_for_y_tol().present()) optimizerIPOPT.alpha_for_y_tol = mu.getAlpha_for_y_tol().get();

	optimizerIPOPT.recalc_y = mu.getRecalc_yDefaultValue();
	if (mu.getRecalc_y().present()) optimizerIPOPT.recalc_y = mu.getRecalc_y().get();

	optimizerIPOPT.recalc_y_feas_tol = mu.getRecalc_y_feas_tolDefaultValue();
	if (mu.getRecalc_y_feas_tol().present()) optimizerIPOPT.recalc_y_feas_tol = mu.getRecalc_y_feas_tol().get();

	// line search
	IPOPTComplexType:: LineSearchType lsc = ipopt.getLineSearch();

	optimizerIPOPT.max_soc = lsc.getMax_socDefaultValue();
	if (lsc.getMax_soc().present()) optimizerIPOPT.max_soc = lsc.getMax_soc().get();

	optimizerIPOPT.watchdog_shortened_iter_trigger = lsc.getWatchdog_shortened_iter_triggerDefaultValue();
	if (lsc.getWatchdog_shortened_iter_trigger().present()) optimizerIPOPT.watchdog_shortened_iter_trigger = lsc.getWatchdog_shortened_iter_trigger().get();

	optimizerIPOPT.watchdog_trial_iter_max = lsc.getWatchdog_trial_iter_maxDefaultValue();
	if (lsc.getWatchdog_trial_iter_max().present()) optimizerIPOPT.watchdog_trial_iter_max = lsc.getWatchdog_trial_iter_max().get();

	optimizerIPOPT.accept_every_trial_step = lsc.getAccept_every_trial_stepDefaultValue();
	if (lsc.getAccept_every_trial_step().present()) optimizerIPOPT.accept_every_trial_step = lsc.getAccept_every_trial_step().get();

	// initialization
	IPOPTComplexType:: InitializationType init = ipopt.getInitialization();

	optimizerIPOPT.bound_frac = init.getBound_fracDefaultValue();
	if (init.getBound_frac().present()) optimizerIPOPT.bound_frac = init.getBound_frac().get();

	optimizerIPOPT.bound_push = init.getBound_pushDefaultValue();
	if (init.getBound_push().present()) optimizerIPOPT.bound_push = init.getBound_push().get();

	optimizerIPOPT.slack_bound_frac = init.getSlack_bound_fracDefaultValue();
	if (init.getSlack_bound_frac().present()) optimizerIPOPT.slack_bound_frac = init.getSlack_bound_frac().get();

	optimizerIPOPT.slack_bound_push = init.getSlack_bound_pushDefaultValue();
	if (init.getSlack_bound_push().present()) optimizerIPOPT.slack_bound_push = init.getSlack_bound_push().get();

	optimizerIPOPT.bound_mult_init_val = init.getBound_mult_init_valDefaultValue();
	if (init.getBound_mult_init_val().present()) optimizerIPOPT.bound_mult_init_val = init.getBound_mult_init_val().get();

	optimizerIPOPT.constr_mult_init_max = init.getConstr_mult_init_maxDefaultValue();
	if (init.getConstr_mult_init_max().present()) optimizerIPOPT.constr_mult_init_max = init.getConstr_mult_init_max().get();

	optimizerIPOPT.bound_mult_init_method = init.getBound_mult_init_methodDefaultValue();
	if (init.getBound_mult_init_method().present()) optimizerIPOPT.bound_mult_init_method = init.getBound_mult_init_method().get();

	// warm start
	IPOPTComplexType:: WarmStartType ws = ipopt.getWarmStart();

	optimizerIPOPT.warm_start_init_point = ws.getWarm_start_init_pointDefaultValue();
	if (ws.getWarm_start_init_point().present()) optimizerIPOPT.warm_start_init_point = ws.getWarm_start_init_point().get();

	optimizerIPOPT.warm_start_bound_frac = ws.getWarm_start_bound_fracDefaultValue();
	if (ws.getWarm_start_bound_frac().present()) optimizerIPOPT.warm_start_bound_frac = ws.getWarm_start_bound_frac().get();

	optimizerIPOPT.warm_start_bound_push = ws.getWarm_start_bound_pushDefaultValue();
	if (ws.getWarm_start_bound_push().present()) optimizerIPOPT.warm_start_bound_push = ws.getWarm_start_bound_push().get();

	optimizerIPOPT.warm_start_slack_bound_frac = ws.getWarm_start_slack_bound_fracDefaultValue();
	if (ws.getWarm_start_slack_bound_frac().present()) optimizerIPOPT.warm_start_slack_bound_frac = ws.getWarm_start_slack_bound_frac().get();

	optimizerIPOPT.warm_start_slack_bound_push = ws.getWarm_start_slack_bound_pushDefaultValue();
	if (ws.getWarm_start_slack_bound_push().present()) optimizerIPOPT.warm_start_slack_bound_push = ws.getWarm_start_slack_bound_push().get();

	optimizerIPOPT.warm_start_mult_bound_push = ws.getWarm_start_mult_bound_pushDefaultValue();
	if (ws.getWarm_start_mult_bound_push().present()) optimizerIPOPT.warm_start_mult_bound_push = ws.getWarm_start_mult_bound_push().get();

	optimizerIPOPT.warm_start_mult_init_max = ws.getWarm_start_mult_init_maxDefaultValue();
	if (ws.getWarm_start_mult_init_max().present()) optimizerIPOPT.warm_start_mult_init_max = ws.getWarm_start_mult_init_max().get();

	// barrier parameter
	IPOPTComplexType::BarrierParameterType bp = ipopt.getBarrierParameter();

	optimizerIPOPT.mehrotra_algorithm = bp.getMehrotra_algorithmDefaultValue();
	if (bp.getMehrotra_algorithm().present()) optimizerIPOPT.mehrotra_algorithm = bp.getMehrotra_algorithm().get();

	optimizerIPOPT.mu_strategy = bp.getMu_strategyDefaultValue();
	if (bp.getMu_strategy().present()) optimizerIPOPT.mu_strategy = bp.getMu_strategy().get();

	optimizerIPOPT.mu_oracle = bp.getMu_oracleDefaultValue();
	if (bp.getMu_oracle().present()) optimizerIPOPT.mu_oracle = bp.getMu_oracle().get();

	optimizerIPOPT.quality_function_max_section_steps = bp.getQuality_function_max_section_stepsDefaultValue();
	if (bp.getQuality_function_max_section_steps().present()) optimizerIPOPT.quality_function_max_section_steps = bp.getQuality_function_max_section_steps().get();

	optimizerIPOPT.fixed_mu_oracle = bp.getFixed_mu_oracleDefaultValue();
	if (bp.getFixed_mu_oracle().present()) optimizerIPOPT.fixed_mu_oracle = bp.getFixed_mu_oracle().get();

	optimizerIPOPT.adaptive_mu_globalization = bp.getAdaptive_mu_globalizationDefaultValue();
	if (bp.getAdaptive_mu_globalization().present()) optimizerIPOPT.adaptive_mu_globalization = bp.getAdaptive_mu_globalization().get();

	optimizerIPOPT.mu_init = parInt->getDblParameter(bp.getMu_initDefaultValue());
	if (bp.getMu_init().present()) optimizerIPOPT.mu_init = parInt->getDblParameter(bp.getMu_init().get());

	optimizerIPOPT.mu_max_fact = bp.getMu_max_factDefaultValue();
	if (bp.getMu_max_fact().present()) optimizerIPOPT.mu_max_fact = bp.getMu_max_fact().get();

	optimizerIPOPT.mu_max = bp.getMu_maxDefaultValue();
	if (bp.getMu_max().present()) optimizerIPOPT.mu_max = bp.getMu_max().get();

	optimizerIPOPT.mu_min = bp.getMu_minDefaultValue();
	if (bp.getMu_min().present()) optimizerIPOPT.mu_min = bp.getMu_min().get();

	optimizerIPOPT.mu_target = bp.getMu_targetDefaultValue();
	if (bp.getMu_target().present()) optimizerIPOPT.mu_target = bp.getMu_target().get();

	optimizerIPOPT.barrier_tol_factor = bp.getBarrier_tol_factorDefaultValue();
	if (bp.getBarrier_tol_factor().present()) optimizerIPOPT.barrier_tol_factor = bp.getBarrier_tol_factor().get();

	optimizerIPOPT.mu_linear_decrease_factor = bp.getMu_linear_decrease_factorDefaultValue();
	if (bp.getMu_linear_decrease_factor().present()) optimizerIPOPT.mu_linear_decrease_factor = bp.getMu_linear_decrease_factor().get();

	optimizerIPOPT.mu_superlinear_decrease_power = bp.getMu_superlinear_decrease_powerDefaultValue();
	if (bp.getMu_superlinear_decrease_power().present()) optimizerIPOPT.mu_superlinear_decrease_power = bp.getMu_superlinear_decrease_power().get();

	// quasi newton
	IPOPTComplexType::QuasiNewtonType qn = ipopt.getQuasiNewton();

	optimizerIPOPT.limited_memory_max_history = parInt->getIntParameter(qn.getLimited_memory_max_historyDefaultValue());
	if (qn.getLimited_memory_max_history().present()) optimizerIPOPT.limited_memory_max_history = parInt->getIntParameter(qn.getLimited_memory_max_history().get());

	optimizerIPOPT.limited_memory_max_skipping = parInt->getIntParameter(qn.getLimited_memory_max_skippingDefaultValue());
	if (qn.getLimited_memory_max_skipping().present()) optimizerIPOPT.limited_memory_max_skipping = parInt->getIntParameter(qn.getLimited_memory_max_skipping().get());


	// derivative checker
	IPOPTComplexType::DerivativeCheckerType dcheck = ipopt.getDerivativeChecker();

	optimizerIPOPT.derivative_test = dcheck.getDerivative_testDefaultValue();
	if (dcheck.getDerivative_test().present()) optimizerIPOPT.derivative_test = dcheck.getDerivative_test().get();

	optimizerIPOPT.derivative_test_perturbation = dcheck.getDerivative_test_perturbationDefaultValue();
	if (dcheck.getDerivative_test_perturbation().present()) optimizerIPOPT.derivative_test_perturbation = dcheck.getDerivative_test_perturbation().get();

	optimizerIPOPT.derivative_test_tol = dcheck.getDerivative_test_tolDefaultValue();
	if (dcheck.getDerivative_test_tol().present()) optimizerIPOPT.derivative_test_tol = dcheck.getDerivative_test_tol().get();

	optimizerIPOPT.derivative_test_print_all = dcheck.getDerivative_test_print_allDefaultValue();
	if (dcheck.getDerivative_test_print_all().present()) optimizerIPOPT.derivative_test_print_all = dcheck.getDerivative_test_print_all().get();

	optimizerIPOPT.jacobian_approximation = dcheck.getJacobian_approximationDefaultValue();
	if (dcheck.getJacobian_approximation().present()) optimizerIPOPT.jacobian_approximation = dcheck.getJacobian_approximation().get();

	optimizerIPOPT.findiff_perturbation = dcheck.getFindiff_perturbationDefaultValue();
	if (dcheck.getFindiff_perturbation().present()) optimizerIPOPT.findiff_perturbation = dcheck.getFindiff_perturbation().get();

	return optimizerIPOPT;
}


rtcRuntimeConfigSettings::SA getSASettings(fews::OptimizerComplexType::SAType sa, parameterInterface *parInt)
{
	rtcRuntimeConfigSettings::SA optimizerSA;

	optimizerSA.ftoll = parInt->getDblParameter(sa.getFtollDefaultValue());
	if (sa.getFtoll().present()) optimizerSA.ftoll = parInt->getDblParameter(sa.getFtoll().get());

	optimizerSA.dels = parInt->getDblParameter(sa.getDelsDefaultValue());
	if (sa.getDels().present()) optimizerSA.dels = parInt->getDblParameter(sa.getDels().get());

	optimizerSA.max_iter = (int)sa.getMax_iterDefaultValue();
	if (sa.getMax_iter().present()) optimizerSA.max_iter = (int)sa.getMax_iter().get();

	optimizerSA.max_cpu_time = parInt->getDblParameter(sa.getMax_cpu_timeDefaultValue());
	if (sa.getMax_cpu_time().present()) optimizerSA.max_cpu_time = parInt->getDblParameter(sa.getMax_cpu_time().get());

	return optimizerSA;
}


rtcRuntimeConfigSettings::SNOPT getSNOPTSettings(fews::OptimizerComplexType::SNOPTType snopt, parameterInterface *parInt)
{
	rtcRuntimeConfigSettings::SNOPT optimizerSNOPT;

	return optimizerSNOPT;
}


void rtcToolsRuntime::parseRuntimeConfigFile()
{
	// xsd file
	::xml_schema::Properties rtcRunTimeConfigProperties;
    rtcRunTimeConfigProperties.schema_location("http://www.wldelft.nl/fews",
		utils::xsd_filename(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "rtcRuntimeConfig.xsd")));

	// xml file
	string filename = utils::getAbsoluteFilename(runtimeSettings.workDir, "rtcRuntimeConfig.xml");
	assertFile(filename);

    // parsing file with runtime information
	unique_ptr<RtcRuntimeConfigComplexType> root;
	try {
		root = parseRtcRuntimeConfig(filename, 0, rtcRunTimeConfigProperties);
	} catch (const xml_schema::Exception &e) {
		cout << e << endl;
		throw runtime_error("error parsing rtcRuntimeConfig.xml file");
	}

	piDiagInterface::addLine(4, "rtcToolsRuntime: runtimeSettings.schemaDir = " + runtimeSettings.schemaDir);
	piDiagInterface::addLine(4, "rtcToolsRuntime: runtimeSettings.workDir = " + runtimeSettings.workDir);

    // files
    if (root->getFiles().present()) {
		RtcRuntimeConfigComplexType::FilesType ft = root->getFiles().get();
		// required
        runtimeSettings.dataConfigFile = ft.getRtcDataConfig().getName();
		runtimeSettings.toolsConfigFile = ft.getRtcToolsConfig().getName();
		// optional
		if (ft.getRtcObjectiveConfig().present()) {
			runtimeSettings.objectiveConfigFile = ft.getRtcObjectiveConfig().get().getName();
		}
		if (ft.getRtcScenarioTreeConfig().present()) {
			runtimeSettings.scenarioTreeConfigFile = ft.getRtcScenarioTreeConfig().get().getName();
		}
		int nPar = (int)ft.getRtcParameterConfig().size();
		if (nPar>0) {
			runtimeSettings.parameterConfigFile.resize(nPar);
			for (int i=0; i<nPar; i++) {
				rtcRuntimeConfigSettings::parameterFile pf;
				pf.filename = ft.getRtcParameterConfig()[i].getName();
				if (ft.getRtcParameterConfig()[i].getType().compare("TREEVECTOR")==0) {
					pf.type = TREEVECTOR;
				} else if (ft.getRtcParameterConfig()[i].getType().compare("PIMODELPARAMETERS")==0) {
					pf.type = PIMODELPARAMETERS;
				}
				pf.prefix = PREFIX_NONE;
				if (ft.getRtcParameterConfig()[i].getPrefix().compare("LOCATIONID")==0) {
					pf.prefix = PREFIX_LOCATIONID;
				}
				runtimeSettings.parameterConfigFile[i] = pf;
			}
		}
    }

	// import of the parameter file if present
	parInt = new parameterInterface();
	for (int i=0; i<(int)runtimeSettings.parameterConfigFile.size(); i++) {
		filename = utils::getAbsoluteFilename(runtimeSettings.workDir, runtimeSettings.parameterConfigFile[i].filename);
		if (utils::fileAvailable(filename)) {
			piDiagInterface::addLine(4, "rtcToolsRuntime: parameter file available", RTCTOOLSRUNTIME_CODE);
			parInt->addFileContent(runtimeSettings.schemaDir, filename,
				runtimeSettings.parameterConfigFile[i].type,
				runtimeSettings.parameterConfigFile[i].prefix);
		}
	}

	// logging
	runtimeSettings.constraintViolationAsError = false;
	runtimeSettings.constraintViolationTol = 1e-4;

	runtimeSettings.outputMetaData = false;

	runtimeSettings.outputObjectiveFunction = false;
	runtimeSettings.outputObjectiveFunctionGradient = false;
	runtimeSettings.outputConstraint = false;
	runtimeSettings.outputJacobian = false;

    if (root->getLogging().present()) {

		RtcRuntimeConfigComplexType::LoggingType log = root->getLogging().get();

		piDiagInterface::setLogLevel(parInt->getIntParameter(log.getLogLevel()));
		if (log.getEventCode().present()) {
			piDiagInterface::setEventCode(parInt->getBoolParameter(log.getEventCode().get()));
		}
		piDiagInterface::setFlush(parInt->getBoolParameter(log.getFlushing()));

		if (log.getConstraintViolationTolerance().present()) {
			runtimeSettings.constraintViolationTol =
				parInt->getDblParameter(log.getConstraintViolationTolerance().get());
		}

		if (log.getReportConstraintViolation().present()) {
			runtimeSettings.constraintViolationAsError = true;
			runtimeSettings.constraintViolationLevel =
				parInt->getIntParameter(log.getReportConstraintViolation().get().getLevel());
		}

		if (log.getOutputMetaData().present() && log.getOutputMetaData().get()) {
			runtimeSettings.outputMetaData = true;
		}

		if (log.getOutputObjectiveFunction().present()) {

			string s = log.getOutputObjectiveFunction().get();
			if (s.compare("VALUE")==0 || s.compare("VALUE+DERIVATIVE")==0) {
				runtimeSettings.outputObjectiveFunction = true;
			}
			if (s.compare("DERIVATIVE")==0 || s.compare("VALUE+DERIVATIVE")==0) {
				runtimeSettings.outputObjectiveFunctionGradient = true;
			}

			s = log.getOutputConstraints().get();
			if (s.compare("VALUE")==0 || s.compare("VALUE+DERIVATIVE")==0) {
				runtimeSettings.outputConstraint = true;
			}
			if (s.compare("DERIVATIVE")==0 || s.compare("VALUE+DERIVATIVE")==0) {
				runtimeSettings.outputJacobian = true;
			}
		}
    }

	// get period from user definition
	if (root->getPeriod().getUserDefined().present()) {

		RtcRuntimeConfigComplexType::PeriodType::UserDefinedType ud = root->getPeriod().getUserDefined().get();

		// start and end time
		t1 = getTime(ud.getStartDate().getDate(), ud.getStartDate().getTime());
		t2 = getTime(ud.getEndDate().getDate(), ud.getEndDate().getTime());

		dt = getTimeStep(ud.getTimeStep().getUnit());
		dt *= parInt->getIntParameter(ud.getTimeStep().getMultiplier());
		dt /= parInt->getIntParameter(ud.getTimeStep().getDivider());

		// optional forecast time, default is the start time t1
		T0 = t1;
		if (ud.getForecastDate().present()) {
			T0 = getTime(ud.getForecastDate().get().getDate(), ud.getForecastDate().get().getTime());
		}
		iForecastTime = (int)((T0-t1)/dt);

		nEnsemble = 1;
		if (ud.getNumberEnsembles().present()) {
			nEnsemble = (int)ud.getNumberEnsembles().get();
		}
	}

	// get period from PI XML time series file
	if (root->getPeriod().getPIInput().present()) {

		RtcRuntimeConfigComplexType::PeriodType::PIInputType pi = root->getPeriod().getPIInput().get();

		// get required info from pi file
		string filename = utils::getAbsoluteFilename(runtimeSettings.workDir, pi.getFile());
		if (!utils::fileAvailable(filename)) {
			piDiagInterface::addLine(1, "void rtcToolsRuntime::parseRunTimeConfigSettings(), file not found - " + filename, RTCTOOLSRUNTIME_CODE);
			throw runtime_error("void rtcToolsRuntime::parseRunTimeConfigSettings() - error - time series import file not found");
		}
		::xml_schema::Properties timeSeriesCollectionComplexTypeProperties;
		timeSeriesCollectionComplexTypeProperties.schema_location("http://www.wldelft.nl/fews/PI",
			utils::xsd_filename(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_timeseries.xsd")));
		unique_ptr<TimeSeriesCollectionComplexType> TimeSeries;
		try {
			TimeSeries = parseTimeSeries(filename, 0, timeSeriesCollectionComplexTypeProperties);
		} catch (const xml_schema::Exception &e) {
			cout << e << endl;
			throw runtime_error("error parsing timeseries_import.xml file");
		}

		// the first time series available
		PI::TimeSeriesComplexType::HeaderType h = TimeSeries->getSeries()[0].getHeader();

		// start and end time
		t1 = getTime(h.getStartDate().getDate(), h.getStartDate().getTime());
		t2 = getTime(h.getEndDate().getDate(), h.getEndDate().getTime());

		dt = getTimeStep(h.getTimeStep().getUnit());
		dt *= h.getTimeStep().getMultiplier();
		dt /= h.getTimeStep().getDivider();

		// optional forecast time, default is the start time t1
		T0 = t1;
		if (h.getForecastDate().present()) {
			T0 = getTime(h.getForecastDate().get().getDate(), h.getForecastDate().get().getTime());
		}
		iForecastTime = (int)((T0-t1)/dt);

		nEnsemble = 1;
		if (pi.getNumberEnsembles().present()) {
		   nEnsemble = (int)pi.getNumberEnsembles().get();
		}
	}

	// get period from PI XML run file
	if (root->getPeriod().getPIRunFile().present()) {

		RtcRuntimeConfigComplexType::PeriodType::PIRunFileType pi = root->getPeriod().getPIRunFile().get();

		// get required info from pi file
		string filename = utils::getAbsoluteFilename(runtimeSettings.workDir, pi.getFile());
		if (!utils::fileAvailable(filename)) {
			piDiagInterface::addLine(1, "void rtcToolsRuntime::parseRunTimeConfigSettings(), file not found - " + filename, RTCTOOLSRUNTIME_CODE);
			throw runtime_error("void rtcToolsRuntime::parseRunTimeConfigSettings() - error - pi run file not found");
		}
		::xml_schema::Properties runComplexTypeProperties;
		runComplexTypeProperties.schema_location("http://www.wldelft.nl/fews/PI",
			utils::xsd_filename(utils::getAbsoluteFilename(runtimeSettings.schemaDir, "pi_run.xsd")));
		unique_ptr<RunComplexType> run;
		try {
			run = parseRun(filename, 0, runComplexTypeProperties);
		} catch (const xml_schema::Exception &e) {
			cout << e << endl;
			throw runtime_error("error parsing run.xml file");
		}

		// start and end time
		t1 = getTime(run->getStartDateTime().getDate(), run->getStartDateTime().getTime());
		t2 = getTime(run->getEndDateTime().getDate(), run->getEndDateTime().getTime());

		dt = getTimeStep(pi.getTimeStep().getUnit());
		dt *= parInt->getIntParameter(pi.getTimeStep().getMultiplier());
		dt /= parInt->getIntParameter(pi.getTimeStep().getDivider());

		// optional forecast time, default is the start time t1
		T0 = getTime(run->getTime0().getDate(), run->getTime0().getTime());
		iForecastTime = (int)((T0-t1)/dt);

		nEnsemble = 1;
		if (pi.getNumberEnsembles().present()) {
		   nEnsemble = (int)pi.getNumberEnsembles().get();
		}
	}

	char buffer[500];
    char b[4][50];
    snprintf(buffer, sizeof(buffer),"%s %s - %s %s, dt[ms] = %d",
		utils::time2datestring(t1, b[0]), utils::time2timestring(t1, b[1]),
		utils::time2datestring(t2, b[2]), utils::time2timestring(t2, b[3]), dt);
    piDiagInterface::addLine(3, "rtcToolsRuntime: simulation period = " + string(buffer), RTCTOOLSRUNTIME_CODE);

    snprintf(buffer, sizeof(buffer), "%d", nEnsemble);
    piDiagInterface::addLine(4, "rtcToolsRuntime: number of ensemble members = " + string(buffer), RTCTOOLSRUNTIME_CODE);

	this->executeObjectiveFunction = true;
	this->executeConstraints = true;

    // single mode definition
    if (root->getMode().present()) {

        runtimeSettings.modeInfo.push_back(getModeInfo(root->getMode().get()));

	// definition of multiple modes
    } else if (root->getModes().present()) {

		for (int i=0; i<(int)root->getModes().get().getMode().size(); i++) {
			runtimeSettings.modeInfo.push_back(
				getModeInfo(root->getModes().get().getMode()[i]));
		}

	// defaut settings
	} else {

		rtcRuntimeConfigSettings::modeInfoContainer c;

		// default is a simulation over the complete period
		c.mode = SIMULATE;
		c.period = COMPLETE;

		runtimeSettings.modeInfo.push_back(c);
	}

	// parallelization
	runtimeSettings.nThread = -1;
	runtimeSettings.parallelEnsembleSim = false;
	runtimeSettings.parallelEnsembleCon = false;
	runtimeSettings.parallelInternalSim = false;
	runtimeSettings.parallelInternalCon = false;

	if (root->getParallelization().present()) {

		// number of maximum threads
		if (root->getParallelization().get().getNThread().present()) {
			runtimeSettings.nThread = parInt->getIntParameter(root->getParallelization().get().getNThread().get());
		}

		// parallelization options for simulation and objective function evaluation
		if (root->getParallelization().get().getSimulation().present()) {
			if (root->getParallelization().get().getSimulation().get().find("ENSEMBLE") != string::npos) {
				runtimeSettings.parallelEnsembleSim = true;
			}
			if (root->getParallelization().get().getSimulation().get().find("INTERNAL") != string::npos) {
				runtimeSettings.parallelInternalSim = true;
			}
		}

		// parallelization options for constraints and constaint Jacobian
		if (root->getParallelization().get().getConstraints().present()) {
			if (root->getParallelization().get().getConstraints().get().find("ENSEMBLE") != string::npos) {
				runtimeSettings.parallelEnsembleCon = true;
			}
			if (root->getParallelization().get().getConstraints().get().find("INTERNAL") != string::npos) {
				runtimeSettings.parallelInternalCon = true;
			}
		}
	}
}


rtcRuntimeConfigSettings::modeInfoContainer rtcToolsRuntime::getModeInfo(fews::RtcRuntimeConfigComplexType::ModeType m) {

	rtcRuntimeConfigSettings::modeInfoContainer c;

    if (m.getSimulation().present()) {

		c.mode = SIMULATE;
		if (m.getSimulation().get().getLimitedMemory().present()) {
			runtimeSettings.limitedMemory = m.getSimulation().get().getLimitedMemory().get();
		}

		c.period = COMPLETE;
		if (m.getSimulation().get().getPeriod().present()) {
			if (m.getSimulation().get().getPeriod().get().compare("UPDATE")==0) {
				c.period = UPDATE;
			}
			if (m.getSimulation().get().getPeriod().get().compare("FORECAST")==0) {
				c.period = FORECAST;
			}
		}

		c.executeObjectiveFunction = true;
		c.executeConstraints = true;
		if (m.getSimulation().get().getExecuteObjectiveFunction().present()) {
			if (!m.getSimulation().get().getExecuteObjectiveFunction().get()) c.executeObjectiveFunction = false;
		}
		if (m.getSimulation().get().getExecuteConstraints().present()) {
			if (!m.getSimulation().get().getExecuteConstraints().get()) c.executeConstraints = false;
		}

	} else if (m.getFirstOrderSensitivity().present()) {

		c.mode = FIRSTORDERSENSITIVITY;
		runtimeSettings.limitedMemory = false;
		c.executeObjectiveFunction = true;
		c.executeConstraints = false;

		c.period = COMPLETE;
		if (m.getFirstOrderSensitivity().get().getPeriod().present()) {
			if (m.getFirstOrderSensitivity().get().getPeriod().get().compare("UPDATE")==0) {
				c.period = UPDATE;
			}
			if (m.getFirstOrderSensitivity().get().getPeriod().get().compare("FORECAST")==0) {
				c.period = FORECAST;
			}
		}

	} else if (m.getFirstOrderDerivativeCheck().present()) {

        c.mode = FIRSTORDERDERIVATIVECHECK;
		runtimeSettings.limitedMemory = false;
		c.executeObjectiveFunction = true;
		c.executeConstraints = true;

		c.period = COMPLETE;
		if (m.getFirstOrderDerivativeCheck().get().getPeriod().present()) {
			if (m.getFirstOrderDerivativeCheck().get().getPeriod().get().compare("UPDATE")==0) {
				c.period = UPDATE;
			}
			if (m.getFirstOrderDerivativeCheck().get().getPeriod().get().compare("FORECAST")==0) {
				c.period = FORECAST;
			}
		}

		c.derivativeCheck_method = CENTRAL;
		if (m.getFirstOrderDerivativeCheck().get().getDerivative_test_method().compare("RIDDERS")==0) {
			c.derivativeCheck_method = RIDDERS;
		}
		c.derivativeCheck_pertubation = m.getFirstOrderDerivativeCheck().get().getDerivative_test_perturbation();
		c.derivativeCheck_tol = m.getFirstOrderDerivativeCheck().get().getDerivative_test_tol();
		c.derivativeCheck_exit_on_error = m.getFirstOrderDerivativeCheck().get().getDerivative_test_exit_on_error();

    } else if (m.getOptimization().present()) {

        c.mode = OPTIMIZE;
		runtimeSettings.limitedMemory = false;
		c.executeObjectiveFunction = true;
		c.executeConstraints = true;

		c.period = COMPLETE;
		if (m.getOptimization().get().getPeriod().present()) {
			if (m.getOptimization().get().getPeriod().get().compare("UPDATE")==0) {
				c.period = UPDATE;
			}
			if (m.getOptimization().get().getPeriod().get().compare("FORECAST")==0) {
				c.period = FORECAST;
			}
		}

		// GAMS
		if (m.getOptimization().get().getOptimizer()[0].getGAMS().present()) {

			c.optimizer = GAMS;
			fews::OptimizerComplexType::GAMSType gams = m.getOptimization().get().getOptimizer()[0].getGAMS().get();
			c.optimizerGAMS = getGAMSSettings(gams, parInt);

		// IPOPT
		} else if (m.getOptimization().get().getOptimizer()[0].getIPOPT().present()) {

			c.optimizer = IPOPT;
			fews::OptimizerComplexType::IPOPTType ipopt = m.getOptimization().get().getOptimizer()[0].getIPOPT().get();
			c.optimizerIPOPT = getIPOPTSettings(ipopt, parInt, info);

		// SA
		} else if (m.getOptimization().get().getOptimizer()[0].getSA().present()) {

			c.optimizer = SA;
			fews::OptimizerComplexType::SAType sa = m.getOptimization().get().getOptimizer()[0].getSA().get();
			c.optimizerSA = getSASettings(sa, parInt);

		// SNOPT
		} else if (m.getOptimization().get().getOptimizer()[0].getSNOPT().present()) {

			c.optimizer = SNOPT;
			fews::OptimizerComplexType::SNOPTType snopt = m.getOptimization().get().getOptimizer()[0].getSNOPT().get();
			c.optimizerSNOPT = getSNOPTSettings(snopt, parInt);
		}

	} else if (m.getClosedLoop().present()) {

        c.mode = CLOSEDLOOP;
		runtimeSettings.limitedMemory = false;
        runtimeSettings.forecastHorizon = (int)m.getClosedLoop()->getForecastHorizon();
        runtimeSettings.recedingHorizon = (int)m.getClosedLoop()->getRecedingHorizon();
		c.executeObjectiveFunction = true;
		c.executeConstraints = true;

		// GAMS
		if (m.getClosedLoop().get().getOptimization().getOptimizer()[0].getGAMS().present()) {

			c.optimizer = GAMS;
			fews::OptimizerComplexType::GAMSType gams = m.getClosedLoop().get().getOptimization().getOptimizer()[0].getGAMS().get();
			c.optimizerGAMS = getGAMSSettings(gams, parInt);

		// IPOPT
		} else if (m.getClosedLoop().get().getOptimization().getOptimizer()[0].getIPOPT().present()) {

			c.optimizer = IPOPT;
			fews::OptimizerComplexType::IPOPTType ipopt = m.getClosedLoop().get().getOptimization().getOptimizer()[0].getIPOPT().get();
			c.optimizerIPOPT = getIPOPTSettings(ipopt, parInt, info);

		// SA
		} else if (m.getClosedLoop().get().getOptimization().getOptimizer()[0].getSA().present()) {

			c.optimizer = SA;
			fews::OptimizerComplexType::SAType sa = m.getClosedLoop().get().getOptimization().getOptimizer()[0].getSA().get();
			c.optimizerSA = getSASettings(sa, parInt);

		// SNOPT
		} else if (m.getClosedLoop().get().getOptimization().getOptimizer()[0].getSNOPT().present()) {

			c.optimizer = SNOPT;
			fews::OptimizerComplexType::SNOPTType snopt = m.getClosedLoop().get().getOptimization().getOptimizer()[0].getSNOPT().get();
			c.optimizerSNOPT = getSNOPTSettings(snopt, parInt);
		}

	} else if (m.getPostprocessing().present()) {

		c.mode = POSTPROCESSING;
		c.executeObjectiveFunction = false;
		c.executeConstraints = false;

	} else {
		throw runtime_error("void rtcToolsRuntime::parseRuntimeConfigFile() - error");
	}

	return c;
}


void rtcToolsRuntime::parseInputFiles()
{
	string sDir = runtimeSettings.schemaDir;
	string wDir = runtimeSettings.workDir;

	// check if one mode is different than simulation, in this case allocate adjoint matrix
	bool allocateAdjoint = false;
	for (int i=0; i<(int)runtimeSettings.modeInfo.size(); i++) {
		if (runtimeSettings.modeInfo[i].mode!=SIMULATE) allocateAdjoint = true;
	}

	// read time series data
	string filename = utils::getAbsoluteFilename(wDir, runtimeSettings.dataConfigFile);
	assertFile(filename);
	tsModel = new timeSeriesModel(sDir, wDir, filename, t1, t2, dt, nEnsemble, ensembleMap, allocateAdjoint, runtimeSettings.limitedMemory);
	int nTimeStep = tsModel->getTimeSeriesTensor()->getNTimeStep();

	// read optional scenario tree
	treeInt = 0;
	ensembleMap = vector<int>(nEnsemble);
	for (int i=0; i<nEnsemble; i++) ensembleMap[i] = i;
	filename = utils::getAbsoluteFilename(wDir, runtimeSettings.scenarioTreeConfigFile);
	if (utils::fileAvailable(filename)) {
		piDiagInterface::addLine(4, "rtcToolsRuntime: scenario tree available", RTCTOOLSRUNTIME_CODE);
		treeInt = new scenarioTreeInterface(sDir, filename, tsModel->getTimeSeriesTensor());
	}

	// read schematisation
	filename = utils::getAbsoluteFilename(wDir, runtimeSettings.toolsConfigFile);
	assertFile(filename);
	schema = new schematisation(sDir, filename, tsModel->getTimeSeriesTensor(), parInt, &runtimeSettings);

	// read optional objective function and allocate a tensor for storing the individual objective function term contributions
	obj = new objectiveFunction();
	filename = utils::getAbsoluteFilename(wDir, runtimeSettings.objectiveConfigFile);
	if (utils::fileAvailable(filename)) {
		piDiagInterface::addLine(4, "rtcToolsRuntime: objective function available", RTCTOOLSRUNTIME_CODE);
		obj = new objectiveFunction(sDir, filename,
			nTimeStep, tsModel->getTimeSeriesTensor(), treeInt, schema, parInt);
	}

	// state import
	tsModel->read();

	iStart = 1;
	iEnd = tsModel->getTimeSeriesTensor()->getNTimeStep()-1;

	//set optimizers
	optimizers = vector<rtcToolsOptimizer*>((int)runtimeSettings.modeInfo.size());
	for (int i=0; i<(int)runtimeSettings.modeInfo.size(); i++) {

		if ((runtimeSettings.modeInfo[i].mode == OPTIMIZE) ||
			(runtimeSettings.modeInfo[i].mode == CLOSEDLOOP)) {

			if (runtimeSettings.modeInfo[i].optimizer == GAMS) {
				optimizers[i] = new rtcToolsGAMS_OPT(this, runtimeSettings.modeInfo[i].optimizerGAMS, runtimeSettings.workDir);
			} else if (runtimeSettings.modeInfo[i].optimizer == IPOPT) {
				optimizers[i] = new rtcToolsIPOPT(this, runtimeSettings.modeInfo[i].optimizerIPOPT, runtimeSettings.workDir);
			} else if (runtimeSettings.modeInfo[i].optimizer == SA) {
				optimizers[i] = new rtcToolsSA(this, runtimeSettings.modeInfo[i].optimizerSA);
			} else if (runtimeSettings.modeInfo[i].optimizer == SNOPT) {
#ifdef HAVE_SNOPT
				optimizers[i] = new rtcToolsSNOPT(this, runtimeSettings.workDir);
#else
				throw runtime_error("void rtcToolsRuntime::initialize() - SNOPT support not compiled in");
#endif
			} else {
				throw runtime_error("void rtcToolsRuntime::initialize() - optimizer not implemented");
			}
		}

		if (runtimeSettings.modeInfo[i].mode == POSTPROCESSING) {
			filename = utils::getAbsoluteFilename(wDir, runtimeSettings.postprocessingConfigFile);
			if (utils::fileAvailable(filename)) {
				piDiagInterface::addLine(4, "rtcToolsRuntime: post processing available", RTCTOOLSRUNTIME_CODE);
				optimizers[i] = new postprocessing(sDir, filename, this);
			}
		}
	}

   	// prepare() is already called in constructor of OpenMIINterface, is this call redundant?
	tsModel->getOpenMIInterface()->prepare();

	openOutput();

	// data validation of the input time series
	tsModel->getTimeSeriesTensor()->validate(iForecastTime, true);
}


/**
  * get and set functions
  **/
void rtcToolsRuntime::setParameters(int n, string* id, double* input)
{
	try {
		// put parameters into interface
		this->parInt->setDblParameters(n, id, input);

		// re-initialize schematization with new parameters
		this->schema->initialize(this->parInt);
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "rtcToolsRuntime::setParameters(int n, double *input) - unknown error.", RTCTOOLSRUNTIME_CODE);
		throw;
	}
}


void rtcToolsRuntime::getInput(int n, double *x)
{
	if (n==0) return;

	timeSeriesTensorInterface* tsTensor = tsModel->getTimeSeriesTensor();

	// number of variables, simulation time steps
	int nStepSim = iEnd-iStart+1;
	vector<double> auxArray(nStepSim);

	for (int i=0; i<(int)ensembleMap.size(); i++) {

		// loop over variables
		for (int j=0; j<(int)obj->getNVariable(); j++) {

			objectiveFunction::variable* v = obj->getVariable(j);

			// extract data from value matrix
			double sFactor = v->scalingFactor;
			for (int k=iStart; k<=iEnd; k++) {
				auxArray[k-iStart] = tsTensor->getValue(ensembleMap[i], k, v->iVariableTSM)*sFactor;
			}

			// aggregation
			aggregation* a = v->vAggregation;
			int nStepIn = nStepSim;
			vector<double> auxArray2 = auxArray;
			if (a!=0) {
				nStepIn = a->getNVariable();
				auxArray2.resize(nStepIn);
				a->aggregate(auxArray, auxArray2, false);
			}

			if (v->ensembleMode == objectiveFunction::JOINT) {
				// joint (default) use of one trajectory for every branch if nEnsemble>1
				for (int k=0; k<nStepIn; k++) {
					x[k + v->nStart] = auxArray2[k];
				}
			} else if (v->ensembleMode == objectiveFunction::TREE) {
				// scenario tree mode
				treeInt->tree2total(i, treeInt->getNStepTotal(), (double*)(x + v->nStart), auxArray2);

			} else if (v->ensembleMode == objectiveFunction::INDEPENDENT) {
				// independent variables on every ensemble branch
				for (int k=0; k<nStepIn; k++) {
					x[k + i*nStepIn + v->nStart] = auxArray2[k];
				}
			}
		}
	}
}


void rtcToolsRuntime::setInput(int n, const double *x)
{
	if (n==0) return;

	timeSeriesTensorInterface* tsTensor = tsModel->getTimeSeriesTensor();

	// number of variables, simulation time steps
	int nStepSim = iEnd-iStart+1;
	vector<double> auxArray(nStepSim);

	for (int i=0; i<(int)ensembleMap.size(); i++) {

		for (int j=0; j<(int)obj->getNVariable(); j++) {

			objectiveFunction::variable* v = obj->getVariable(j);
			int index = v->iVariableTSM;

			// aggregation
			aggregation* a = v->vAggregation;
			int nStepIn = nStepSim;
			vector<double> auxArray2 = auxArray;
			if (a!=0) {
				nStepIn = a->getNVariable();
				auxArray2.resize(nStepIn);
			}

			if (v->ensembleMode == objectiveFunction::JOINT) {
				// joint (default) use of one trajectory for every branch if nEnsemble>1
				for (int k=0; k<nStepIn; k++) {
					auxArray2[k] = x[k + v->nStart];
				}
			} else if (v->ensembleMode == objectiveFunction::TREE) {
				// scenario tree mode
				treeInt->total2tree(i, treeInt->getNStepTotal(), (double*)(x + v->nStart), auxArray2);
			} else if (v->ensembleMode == objectiveFunction::INDEPENDENT) {
				// independent variables on every ensemble branch
				for (int k=0; k<nStepIn; k++) {
					auxArray2[k] = x[k + i*nStepIn + v->nStart];
				}
			}

			if (a!=0) {
				for (int k=0; k<nStepSim; k++) auxArray[k] = std::numeric_limits<double>::quiet_NaN();
				a->deaggregate(auxArray2, tsTensor->getValue(ensembleMap[i], iStart-1, index), auxArray);
			} else {
				auxArray = auxArray2;
			}

			// apply scaling factor and put data into time series tensor
			for (int k=iStart; k<=iEnd; k++) {
				tsTensor->setValue(ensembleMap[i], k, index, auxArray[k-iStart]/v->scalingFactor);
			}
		}
	}
}


int rtcToolsRuntime::getN()
{
	int n = 0;

	for (int i=0; i<obj->getNVariable(); i++) {
		n += obj->getVariable(i)->n;
	}

	return n;
}


void rtcToolsRuntime::getN_metadata(vector<metadata> &var_md)
{
	for (int i=0; i<obj->getNVariable(); i++) {
		objectiveFunction::variable* v = obj->getVariable(i);
		for (int j=0; j<v->n; j++) {
			stringstream ss;
			ss << v->id << "@" << j;
			var_md[v->nStart + j].id = v->id;
			var_md[v->nStart + j].label = ss.str();
		}
	}
}


int rtcToolsRuntime::getM()
{
	int m = 0;

	for (int i=0; i<obj->getNVariableRateOfChange(); i++) {
		m += obj->getVariableRateOfChange(i)->nc;
	}
	for (int i=0; i<obj->getNVariableAverage(); i++) {
		m += obj->getVariableAverage(i)->nc;
	}
	for (int i=0; i<obj->getNVariableEqual(); i++) {
		m += obj->getVariableEqual(i)->nc;
	}
	for (int i=0; i<obj->getNState(); i++) {
		m += obj->getState(i)->nc;
	}
	for (int i=0; i<obj->getNVariableChance(); i++) {
		m += obj->getVariableChance(i)->nc;
	}

	return m;
}


void rtcToolsRuntime::getM_metadata(vector<metadata> &con_md)
{
	for (int i=0; i<obj->getNVariableRateOfChange(); i++) {

		string id = obj->getVariableRateOfChange(i)->id;

		for (int j=0; j<(int)obj->getVariableRateOfChange(i)->rof2DArray.size(); j++) {
			for (int k=0; k<(int)obj->getVariableRateOfChange(i)->rof2DArray[j].size(); k++) {

				rateOfChangeConstraint& c = obj->getVariableRateOfChange(i)->rof2DArray[j][k];
				if (c.iStep.size() == 0)
					continue;

				// id and step
				con_md[c.ncIndx].id = id;
				con_md[c.ncIndx].step = c.iStep[1];

				// label
				stringstream ss;
				ss << id << "@" << c.iStep[1] << "(-" << c.iStep[0] << ")";
				con_md[c.ncIndx].label = ss.str();
			}
		}
	}

	for (int i=0; i<obj->getNVariableAverage(); i++) {

		string id = obj->getVariableAverage(i)->id;

		for (int j=0; j<(int)obj->getVariableAverage(i)->av2DArray.size(); j++) {
			for (int k=0; k<(int)obj->getVariableAverage(i)->av2DArray[j].size(); k++) {

				averageConstraint& c = obj->getVariableAverage(i)->av2DArray[j][k];
				if (c.iStep.size() == 0)
					continue;

				// id and step
				con_md[c.ncIndx].id = id;
				con_md[c.ncIndx].step = c.iStep.back();

				// label
				stringstream ss;
				ss << id << "@" << c.iStep.back() << "(-" << c.iStep.front() << ")";
				con_md[c.ncIndx].label = ss.str();
			}
		}
	}

	for (int i=0; i<obj->getNState(); i++) {

		string id = obj->getState(i)->id;

		for (int j=0; j<(int)obj->getState(i)->state2DArray.size(); j++) {
			for (int k=0; k<(int)obj->getState(i)->state2DArray[j].size(); k++) {

				stateConstraint& c = obj->getState(i)->state2DArray[j][k];

				// id and step
				con_md[c.ncIndx].id = id;
				con_md[c.ncIndx].step = c.iStepResiduum;

				// label
				stringstream ss;
				ss << id << "@" << c.iStepResiduum;
				con_md[c.ncIndx].label = ss.str();
			}
		}
	}
}


int rtcToolsRuntime::getNNZ()
{
	int nnz = 0;

	for (int i=0; i<obj->getNVariableRateOfChange(); i++) {
		nnz += obj->getVariableRateOfChange(i)->nnz;
	}
	for (int i=0; i<obj->getNVariableAverage(); i++) {
		nnz += obj->getVariableAverage(i)->nnz;
	}
	for (int i=0; i<obj->getNVariableEqual(); i++) {
		nnz += obj->getVariableEqual(i)->nnz;
	}
	for (int i=0; i<obj->getNState(); i++) {
		nnz += obj->getState(i)->nnz;
	}
	for (int i=0; i<obj->getNVariableChance(); i++) {
		nnz += obj->getVariableChance(i)->nnz;
	}

	return nnz;
}


void rtcToolsRuntime::get_bounds_info(int n, double* x_l, double* x_u,
				                      int m, double* g_l, double* g_u)
{
	// variable
	for (int i=0; i<obj->getNVariable(); i++) {
		objectiveFunction::variable* var = obj->getVariable(i);
		for (int j=0; j<var->n; j++) {
			assert(var->nStart+j < n);
			x_l[var->nStart+j] = var->x_l[j];
			x_u[var->nStart+j] = var->x_u[j];
		}
	}

	// rate of change constraint
	for (int i=0; i<obj->getNVariableRateOfChange(); i++) {
		objectiveFunction::variableRateOfChange* roc = obj->getVariableRateOfChange(i);
		for (int j=0; j<(int)roc->rof2DArray.size(); j++) {
			for (int k=0; k<(int)roc->rof2DArray[j].size(); k++) {
				assert(roc->rof2DArray[j][k].ncIndx < m);
				g_l[roc->rof2DArray[j][k].ncIndx] = roc->rof2DArray[j][k].g_l;
				g_u[roc->rof2DArray[j][k].ncIndx] = roc->rof2DArray[j][k].g_u;
			}
		}
	}

	// average constraint
	for (int i=0; i<obj->getNVariableAverage(); i++) {
		objectiveFunction::variableAverage* av = obj->getVariableAverage(i);
		for (int j=0; j<(int)av->av2DArray.size(); j++) {
			for (int k=0; k<(int)av->av2DArray[j].size(); k++) {
				assert(av->av2DArray[j][k].ncIndx < m);
				g_l[av->av2DArray[j][k].ncIndx] = av->av2DArray[j][k].g_l;
				g_u[av->av2DArray[j][k].ncIndx] = av->av2DArray[j][k].g_u;
			}
		}
	}

	// equality constraint
	for (int i=0; i<obj->getNVariableEqual(); i++) {
		objectiveFunction::variableEqual* eq = obj->getVariableEqual(i);
		for (int j=0; j<(int)eq->eq2DArray.size(); j++) {
			for (int k=0; k<(int)eq->eq2DArray[j].size(); k++) {
				assert(eq->eq2DArray[j][k].ncIndx < m);
				g_l[eq->eq2DArray[j][k].ncIndx] = 0.0;
				g_u[eq->eq2DArray[j][k].ncIndx] = 0.0;
			}
		}
	}

	// state constraint
	for (int i=0; i<obj->getNState(); i++) {
		objectiveFunction::state* s = obj->getState(i);
		for (int j=0; j<(int)s->state2DArray.size(); j++) {
			for (int k=0; k<(int)s->state2DArray[j].size(); k++) {
				assert(s->state2DArray[j][k].ncIndx < m);
				g_l[s->state2DArray[j][k].ncIndx] = s->scalingFactor * s->state2DArray[j][k].g_l;
				g_u[s->state2DArray[j][k].ncIndx] = s->scalingFactor * s->state2DArray[j][k].g_u;
			}
		}
	}

	// chance constraint
	for (int i=0; i<obj->getNVariableChance(); i++) {
		objectiveFunction::variableChance* cc = obj->getVariableChance(i);
		for (int j=0; j<(int)cc->chance1DArray.size(); j++) {
			assert(cc->chance1DArray[j].ncIndx < m);
			g_l[cc->chance1DArray[j].ncIndx] = cc->chance1DArray[j].g_l;
			g_u[cc->chance1DArray[j].ncIndx] = cc->chance1DArray[j].g_u;
		}
	}

	// some checks
	for (int i=0; i<n; i++) {
		if (x_l[i]==x_l[i] && x_u[i]==x_u[i] && x_u[i]<x_l[i]) {
			// upper bound should not be < lower bound, if both values are set
			vector<metadata> var_md(n);
			getN_metadata(var_md);
			stringstream ss;
			ss << "void rtcToolsRuntime::get_bounds_info(...) - lower bound > upper bound for variable " << var_md[i].label << ". "
			   << "Lower bound: " << setprecision(12) << x_l[i] << ". Upper bound: " << setprecision(12) << x_u[i];
			piDiagInterface::addLine(1, ss.str(), RTCTOOLSRUNTIME_CODE);
			throw runtime_error(ss.str().c_str());
		}
	}

	for (int i=0; i<m; i++) {
		if (g_l[i]==g_l[i] && g_u[i]==g_u[i] && g_u[i]<g_l[i]) {
			// upper bound should not be < lower bound, if both values are set
			vector<metadata> con_md(n);
			getM_metadata(con_md);
			stringstream ss;
			ss << "void rtcToolsRuntime::get_bounds_info(...) - lower bound > upper bound for constraint " << con_md[i].label << ". "
			   << "Lower bound: " << setprecision(12) << g_l[i] << ". Upper bound: " << setprecision(12) << g_u[i];
			piDiagInterface::addLine(1, ss.str(), RTCTOOLSRUNTIME_CODE);
			throw runtime_error(ss.str().c_str());
		}
	}
}


void rtcToolsRuntime::eval_g(int n, const double* x, int m, double* g)
{
	// acting on every ensembe member or scenario tree branch
	if (runtimeSettings.parallelEnsembleCon) {
		#pragma omp parallel for schedule(dynamic)
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].eval_g(n, x, m, g);
		}
	} else {
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].eval_g(n, x, m, g);
		}
	}

	// acting on the COMPLETE ensembe or scenario tree
	// TODO
	double*** ten = tsModel->getTimeSeriesTensor()->getValueTensor();

	// chance constraints
	for (int i=0; i<(int)obj->getNVariableChance(); i++) {
		//
	}

	// check g for valid numbers
	for (int i=0; i<m; i++) {
		if (g[i]!=g[i]) {

			vector<metadata> con_md(m);
			getM_metadata(con_md);

			piDiagInterface::addLine(1, "void rtcToolsRuntime::eval_g() - error - NaN detected in constraint vector at index " + con_md[i].label, RTCTOOLSRUNTIME_CODE);
		}
	}
}


void rtcToolsRuntime::eval_jac_g(int m, int nnz, int* iRow, int* jCol, double* values)
{
	// acting on every ensembe member or scenario tree branch
	if (runtimeSettings.parallelEnsembleCon) {
		#pragma omp parallel for schedule(dynamic)
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].eval_jac_g(m, nnz, iRow, jCol, values);
		}
	} else {
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].eval_jac_g(m, nnz, iRow, jCol, values);
		}
	}

	// acting on the COMPLETE ensembe or scenario tree
	// TODO
	double*** ten = tsModel->getTimeSeriesTensor()->getValueTensor();

	// chance constraints
	for (int i=0; i<(int)obj->getNVariableChance(); i++) {
		//
	}

	// check jac_g for valid numbers
	if (values) {
		for (int i=0; i<nnz; i++) {
			if (values[i]!=values[i]) {

				vector<metadata> var_md(this->getN());
				getN_metadata(var_md);

				vector<metadata> con_md(m);
				getM_metadata(con_md);

				if (iRow != NULL && jCol != NULL)
					piDiagInterface::addLine(1, "void rtcToolsRuntime::eval_g() - error - NaN detected in constraint Jacobian at index " + con_md[iRow[i]].label + "," + var_md[jCol[i]].label, RTCTOOLSRUNTIME_CODE);
				else
					piDiagInterface::addLine(1, "void rtcToolsRuntime::eval_g() - error - NaN detected in constraint Jacobian", RTCTOOLSRUNTIME_CODE);
			}
		}
	}
}


double rtcToolsRuntime::simulate()
{
	try {
		return simulate(0, (double*)0);
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "rtcToolsRuntime::simulate() - error - " + string(e.what()) + ".", RTCTOOLSRUNTIME_CODE);
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "rtcToolsRuntime::simulate() - unknown error.", RTCTOOLSRUNTIME_CODE);
		throw;
	}
}


void rtcToolsRuntime::simulate(int iStep)
{
	try {
		if ((iStep<1) || (iStep>=tsModel->getTimeSeriesTensor()->getNTimeStep())) {
			piDiagInterface::addLine(1, "void rtcToolsRuntime::simulate(int iStep), time step index outside of range.", RTCTOOLSRUNTIME_CODE);
			return;
		}

		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].simulate(iStep);
		}
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "rtcToolsRuntime::simulate(int iStep) - error - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "rtcToolsRuntime::simulate(int iStep) - unknown error in rtcTools", RTCTOOLSRUNTIME_CODE);
		throw;
	}
}


double rtcToolsRuntime::simulate(int n, const double *input, int iTerm, double*** JInc3DArray)
{
	double J = 0;

	vector<double> Je(tsModel->getTimeSeriesTensor()->getEnsembleMap().size());

	if (n>0) setInput(n, input);

	// simulate ensembles and execute the branch-related objective function terms
	if (runtimeSettings.parallelEnsembleSim) {
		#pragma omp parallel for schedule(dynamic)
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			// exception is handled locally to avoid an uncontrolled program termination
			try {
				if (JInc3DArray) {
					Je[i] = rtcSim[i].simulate(iStart, iEnd, this->executeObjectiveFunction, iTerm, JInc3DArray[ensembleMap[i]]);
				} else {
					Je[i] = rtcSim[i].simulate(iStart, iEnd, this->executeObjectiveFunction, iTerm);
				}
			} catch (exception &e) {
				{
					piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - error in parallel simulation - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
					throw runtime_error("double rtcToolsRuntime::simulate(int n, double *input) - error in parallel simulation - " + string(e.what()));
				}
			} catch (...) {
				{
					piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - unknown error in parallel simulation", RTCTOOLSRUNTIME_CODE);
					throw runtime_error("double rtcToolsRuntime::simulate(int n, double *input) - unknown error in parallel simulation - ");
				}
			}
		}
	} else {
		try {
			for (int i=0; i<(int)ensembleMap.size(); i++) {
				if (JInc3DArray) {
					Je[i] = rtcSim[i].simulate(iStart, iEnd, this->executeObjectiveFunction, iTerm, JInc3DArray[ensembleMap[i]]);
				} else {
					Je[i] = rtcSim[i].simulate(iStart, iEnd, this->executeObjectiveFunction, iTerm);
				}
			}
		} catch (exception &e) {
			piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - error in sequential simulation - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
			throw runtime_error("double rtcToolsRuntime::simulate(int n, double *input) - error in sequential simulation - " + string(e.what()));
		} catch (...) {
			piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - unknown error in sequential simulation", RTCTOOLSRUNTIME_CODE);
			throw runtime_error("double rtcToolsRuntime::simulate(int n, double *input) - unknown error in sequential simulation");
		}
	}

	// and add total objective function value
	for (int i=0; i<(int)ensembleMap.size(); i++) {
		J += Je[i];
	}

	if (runtimeSettings.limitedMemory) {
		piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - ensemble-related objective function terms skipped", RTCTOOLSRUNTIME_CODE);
	} else {
		// execute the ensemble-related objective function terms
		for (int i=0; i<(int)obj->getEnsembleTermArray().size(); i++) {
			J += obj->getEnsembleTermArray()[i]->evaluate(iStart, iEnd,
				tsModel->getTimeSeriesTensor()->getValueTensor(), JInc3DArray, ensembleMap, pVec);
		}
	}

	// check for NaN as objective function value
	if (J!=J) {
		piDiagInterface::addLine(1, "double rtcToolsRuntime::simulate(int n, double *input) - objective function value is NaN", RTCTOOLSRUNTIME_CODE);
		throw runtime_error("double rtcToolsRuntime::simulate(int n, double *input) - objective function value is NaN");
	}

	return J;
}


void rtcToolsRuntime::eval_grad_f(int n, double *grad_f, int iTerm)
{
	timeSeriesTensorInterface* tsTensor = tsModel->getTimeSeriesTensor();

	try {
		// initialize adjoint matrices
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			timeseries::timeSeriesMatrix* tsMatrix= dynamic_cast<timeSeriesMatrix*>(tsModel->getTimeSeriesTensor()->getTimeSeriesMatrix(ensembleMap[i]));
			tsMatrix->initializeObj(iStart-1, iEnd, 0.0);
		}

		// execute the ensemble-related objective function terms
		for (int i=0; i<(int)obj->getEnsembleTermArray().size(); i++) {
			obj->getEnsembleTermArray()[i]->evaluateDer(iStart, iEnd, tsTensor->getValueTensor(),
				tsTensor->getObjTensor(), ensembleMap, pVec);
		}

		// execute the branch-related objective function terms and run the adjoint model for each branch
		for (int i=0; i<(int)ensembleMap.size(); i++) {
			rtcSim[i].evaluateGradient(iStart, iEnd, iTerm);
		}

		// initialize gradient vector
		for (int i=0; i<n; i++) grad_f[i] = 0.0;

		// add objective function gradient if required
		if (n==0) return;

		// number of variables, simulation time steps
		int nStepSim = iEnd-iStart+1;
		vector<double> auxArray(nStepSim);

		for (int i=0; i<(int)ensembleMap.size(); i++) {

			//rtcSim[i].getGradient(n, grad_f, iStart, iEnd);
			for (int j=0; j<(int)obj->getNVariable(); j++) {

				objectiveFunction::variable* v = obj->getVariable(j);

				// extract gradient from objective matrix
				for (int k=iStart; k<=iEnd; k++) {
					auxArray[k-iStart] = tsTensor->getObjTensor()[ensembleMap[i]][k][v->iVariableTSM]/v->scalingFactor;
				}

				// aggregation
				aggregation* a = v->vAggregation;
				int nStepIn = nStepSim;
				vector<double> auxArray2 = auxArray;
				if (a!=0) {
					nStepIn = a->getNVariable();
					auxArray2.resize(nStepIn);
					a->aggregate(auxArray, auxArray2, true);
				}

				if (v->ensembleMode == objectiveFunction::JOINT) {
					// joint (default) use of one trajectory for every branch if nEnsemble>1
					for (int k=0; k<nStepIn; k++) {
						grad_f[k + v->nStart] += auxArray2[k];
					}
				} else if (v->ensembleMode == objectiveFunction::TREE) {
					// scenario tree mode
					treeInt->addBranchData(i, treeInt->getNStepTotal(), (double*)(grad_f + v->nStart), auxArray2);

				} else if (v->ensembleMode == objectiveFunction::INDEPENDENT) {
					// independent variables on every ensemble branch
					for (int k=0; k<nStepIn; k++) {
						grad_f[k + i*nStepIn + v->nStart] += auxArray2[k];
					}
				}
			}
		}
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::evaluateGradient(int n, double *dJ) - error - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::evaluateGradient(int n, double *dJ) - unknown error in rtcTools", RTCTOOLSRUNTIME_CODE);
		throw;
	}

	// check grad_f for valid numbers
	for (int i=0; i<n; i++) {
		if (grad_f[i]!=grad_f[i]) {

			vector<metadata> var_md(n);
			getN_metadata(var_md);

			piDiagInterface::addLine(1, "void rtcToolsRuntime::eval_g() - error - NaN detected in gradient vector at index " + var_md[i].label, RTCTOOLSRUNTIME_CODE);
		}
	}
}


double rtcToolsRuntime::eval_grad_f_num(int n, double *x, double EPS, int iX, int iTerm)
{
	const int DIM = 20;
	const double CON = 1.5;
	const double CON2 = CON*CON;
	const double SAFE = 2.0;

	double a[DIM][DIM];
	double EPS2 = EPS;
	double ans;

	double temp = x[iX];

	// central difference
	x[iX] = temp + EPS2;
	double J1 = simulate(n, x, iTerm);
	x[iX] = temp - EPS2;
	double J2 = simulate(n, x, iTerm);
	a[0][0] = (J1-J2) / (2.0*EPS2);

	double err = 1e30;

	for (int i=1; i<DIM; i++) {

		EPS2 /= CON;

		x[iX] = temp + EPS2;
		double J1 = simulate(n, x, iTerm);
		x[iX] = temp - EPS2;
		double J2 = simulate(n, x, iTerm);
		a[0][i] = (J1-J2) / (2.0*EPS2);

		double fac = CON2;

		for (int j=1; j<=i; j++) {
			a[j][i] = (a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
			fac = CON2*fac;
			double errt = max(fabs(a[j][i]-a[j-1][i]), fabs(a[j][i])-a[j-1][i-1]);
			if (errt<=err) {
				err = errt;
				ans = a[j][i];
			}
		}

		if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*err) break;
	}

	x[iX] = temp;

	return ans;
}


bool rtcToolsRuntime::check_grad_f_ridders(int n, double *x, double EPS, double TOL)
{
	vector<metadata> var_md(n);
	getN_metadata(var_md);

	printf("First-order derivative check (Ridders):\n");

	// objective function value and analytical gradient
	double J = simulate(n, x);
	double* grad_f = new double[n];
	eval_grad_f(n, grad_f);

	// numerical gradient
	bool have_errors = false;
	for (int i=0; i<n; i++) {

		//double gn = eval_grad_f_num(n, x, max(EPS, fabs(x[i])*EPS), i);
		double gn = eval_grad_f_num(n, x, EPS, i);
		double diff = fabs(gn-grad_f[i]) / max(1.0, fabs(grad_f[i]));

		if (diff>TOL) {
			printf("grad[%d_%s]: diff= %.3e, num= %.7f, an= %.7f\n", i, var_md[i].label.c_str(), diff, gn, grad_f[i]);

			double* grad_fn = new double[n];

			// further analysis, loop over all objective function terms
			for (int j=0; j<(int)obj->getNTerm(); j++) {

				J = simulate(n, x);
				eval_grad_f(n, grad_fn, j);
				gn = eval_grad_f_num(n, x, EPS, i, j);
				double diff = fabs(gn-grad_fn[i]) / max(1.0, fabs(grad_fn[i]));

				if (diff>0.0) printf(" grad[%s]: diff= %.3e (%.7f, %.7f)\n", obj->getTerm(j)->getID().c_str(), diff, gn, grad_fn[i]);
			}

			delete[] grad_fn;

			have_errors = true;
		}
	}

	printf("\nFinalization of first-order derivative check\n\n\n");

	delete[] grad_f;

	return have_errors;
}


bool rtcToolsRuntime::check_grad_f_central(int n, double *x, double EPS, double TOL)
{
	vector<metadata> var_md(n);
	getN_metadata(var_md);

	printf("First-order derivative check (central):\n");

	// objective function value and analytical gradient
	double J = simulate(n, x);
	double* grad_f = new double[n];
	eval_grad_f(n, grad_f);

	// numerical gradient
	bool have_errors = false;
	for (int i=0; i<n; i++) {

		double temp = x[i];

		// central difference
		double x1 = x[i] = temp - max(fabs(temp)*EPS, EPS);
		double J1 = simulate(n, x);
		double x2 = x[i] = temp + max(fabs(temp)*EPS, EPS);
		double J2 = simulate(n, x);

		double gn = (J2-J1) / (x2-x1);
		double diff = fabs(gn-grad_f[i]) / max(1.0, fabs(grad_f[i]));

		x[i] = temp;

		if (diff>TOL) {
			printf("grad[%d_%s]: diff= %.3e, num= %.7f, an= %.7f\n", i, var_md[i].label.c_str(), diff, gn, grad_f[i]);

			double* grad_fn = new double[n];

			// further analysis, loop over all objective function terms
			for (int j=0; j<(int)obj->getNTerm(); j++) {

				simulate(n, x);
				eval_grad_f(n, grad_fn, j);

				temp = x[i];

				// central difference
				double x1n = x[i] = temp - max(fabs(temp)*EPS, EPS);
				double J1n = simulate(n, x, j);
				double x2n = x[i] = temp + max(fabs(temp)*EPS, EPS);
				double J2n = simulate(n, x, j);

				double gn = (J2n-J1n) / (x2n-x1n);
				diff = fabs(gn-grad_fn[i]) / max(1.0, fabs(grad_fn[i]));
				if (diff>0.0) printf(" grad[%s]: diff= %.3e (%.7f, %.7f)\n", obj->getTerm(j)->getID().c_str(), diff, gn, grad_fn[i]);

				x[i] = temp;
			}

			delete[] grad_fn;

			have_errors = true;
		}
	}

	printf("\nFinalization of first-order derivative check\n\n\n");

	delete[] grad_f;

	return have_errors;
}


void rtcToolsRuntime::openOutput()
{
	try {
		tsModel->openFiles(false);
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::openOutput(void) - error - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::openOutput(void) - unknown error in rtcTools", RTCTOOLSRUNTIME_CODE);
		throw;
	}

	piDiagInterface::write();
}


void rtcToolsRuntime::writeOutput(int timeStep, bool isFinalTimeStep)
{
	try {
		tsModel->write(false, timeStep, isFinalTimeStep);
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::writeOutput(int, bool) - error - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::writeOutput(int, bool) - unknown error in rtcTools", RTCTOOLSRUNTIME_CODE);
		throw;
	}
}


void rtcToolsRuntime::reportMetaData()
{
	int n = getN();
	vector<metadata> var_md(n);
	getN_metadata(var_md);

	// constraints and bounds
	int m = getM();
	vector<metadata> con_md(m);
	getM_metadata(con_md);

	string filename = utils::getAbsoluteFilename(runtimeSettings.workDir, "metadata.txt");
	ofstream txtFile;
	txtFile.open(filename.c_str(), ios::out | ios::trunc);

	for (int i=0; i<n; i++) {
		txtFile << "variable(" << i << ")= " << var_md[i].label << endl;
	}

	for (int i=0; i<m; i++) {
		txtFile << "constraint(" << i << ")= " << con_md[i].label << endl;
	}

	txtFile.close();
}


void rtcToolsRuntime::reportObjectives()
{
	int nObjSeries = 5*obj->getNTerm() + 6*obj->getNEnsembleTerm();
	vector<string> objSeriesCSV = vector<string>(nObjSeries);
	vector<piTimeSeries> objSeriesPI = vector<piTimeSeries>(nObjSeries);

	// meta data for csv and pi xml output
	int indx = 0;
	for (int i=0; i<obj->getNTerm(); i++) {

		objectiveFunctionTerm* o = obj->getTerm(i);

		// csv output
		objSeriesCSV[indx] = o->getID() + "_setpoint";
		objSeriesCSV[indx+1] = o->getID() + "_input";
		objSeriesCSV[indx+2] = o->getID() + "_distance";
		objSeriesCSV[indx+3] = o->getID() + "_j_unweighted";
		objSeriesCSV[indx+4] = o->getID() + "_j_weighted";

		// pi xml output
		string piid = tsModel->getTimeSeriesTensor()->getSeriesIDs()[o->getIInput()];
		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(piid);

		// use the mapping in the rtcDataConfig if available
		// otherwise, fields for location and parameter ID will be empty
		if (pits) {

			// if a pi time series definition is available,
			// copy it, add additional qualifier
			objSeriesPI[indx] = *pits;
			objSeriesPI[indx].addQualifierID("setpoint");
			objSeriesPI[indx].addQualifierID( o->getID());

			objSeriesPI[indx+1] = *pits;
			objSeriesPI[indx+1].addQualifierID("input");
			objSeriesPI[indx+1].addQualifierID( o->getID());

			objSeriesPI[indx+2] = *pits;
			objSeriesPI[indx+2].addQualifierID("distance");
			objSeriesPI[indx+2].addQualifierID( o->getID());

			objSeriesPI[indx+3] = *pits;
			objSeriesPI[indx+3].addQualifierID("j_unweighted");
			objSeriesPI[indx+3].addQualifierID( o->getID());

			objSeriesPI[indx+4] = *pits;
			objSeriesPI[indx+4].addQualifierID("j_weighted");
			objSeriesPI[indx+4].addQualifierID( o->getID());
		}

		// index is set in any case
		objSeriesPI[indx].setIndex(indx);
		objSeriesPI[indx+1].setIndex(indx+1);
		objSeriesPI[indx+2].setIndex(indx+2);
		objSeriesPI[indx+3].setIndex(indx+3);
		objSeriesPI[indx+4].setIndex(indx+4);

		indx += 5;
	}

	for (int i=0; i<obj->getNEnsembleTerm(); i++) {

		ensembleObjectiveFunctionTerm* o = obj->getEnsembleTerm(i);

		// csv output
		objSeriesCSV[indx] = o->getID() + "_bound";
		objSeriesCSV[indx+1] = o->getID() + "_inputMean";
		objSeriesCSV[indx+2] = o->getID() + "_inputStd";
		objSeriesCSV[indx+3] = o->getID() + "_distance";
		objSeriesCSV[indx+4] = o->getID() + "_j_unweighted";
		objSeriesCSV[indx+5] = o->getID() + "_j_weighted";

		// pi xml output
		string piid = tsModel->getTimeSeriesTensor()->getSeriesIDs()[o->getIInput()];
		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(piid);

		// use the mapping in the rtcDataConfig if available
		// otherwise, fields for location and parameter ID will be empty
		if (pits) {

			// if a pi time series definition is available,
			// copy it, add additional qualifier
			objSeriesPI[indx] = *pits;
			objSeriesPI[indx].addQualifierID("bound");
			objSeriesPI[indx].addQualifierID( o->getID());

			objSeriesPI[indx+1] = *pits;
			objSeriesPI[indx+1].addQualifierID("inputMean");
			objSeriesPI[indx+1].addQualifierID( o->getID());

			objSeriesPI[indx+2] = *pits;
			objSeriesPI[indx+2].addQualifierID("inputStd");
			objSeriesPI[indx+2].addQualifierID( o->getID());

			objSeriesPI[indx+3] = *pits;
			objSeriesPI[indx+3].addQualifierID("distance");
			objSeriesPI[indx+3].addQualifierID( o->getID());

			objSeriesPI[indx+4] = *pits;
			objSeriesPI[indx+4].addQualifierID("j_unweighted");
			objSeriesPI[indx+4].addQualifierID( o->getID());

			objSeriesPI[indx+5] = *pits;
			objSeriesPI[indx+5].addQualifierID("j_weighted");
			objSeriesPI[indx+5].addQualifierID( o->getID());
		}

		// index is set in any case
		objSeriesPI[indx].setIndex(indx);
		objSeriesPI[indx+1].setIndex(indx+1);
		objSeriesPI[indx+2].setIndex(indx+2);
		objSeriesPI[indx+3].setIndex(indx+3);
		objSeriesPI[indx+4].setIndex(indx+4);
		objSeriesPI[indx+5].setIndex(indx+5);

		indx += 6;
	}

	// memory allocation for objective function terms (ensembles x time steps x series)
	int nTimeStep = tsModel->getTimeSeriesTensor()->getNTimeStep();
	double*** JInc3D = utils::dten(nEnsemble, nTimeStep, nObjSeries);
	for (int i=0; i<nEnsemble; i++) {
		for (int j=0; j<nTimeStep; j++) {
			for (int k=0; k<nObjSeries; k++) {
				JInc3D[i][j][k] = 0.0;
			}
		}
	}

	// execute objective function with optional array for the objective function terms
	// (for performance reasons, the standard execution is without)
	int n = getN();
	double* x = new double[n];
	getInput(n, x);
	simulate(n, x, -1, JInc3D);

	// csv and pi xml write actions
	if (tsModel->getCsvInterface()) {
		tsModel->getCsvInterface()->writeFiles("objective", objSeriesCSV, JInc3D);
	}
	tsModel->getPISAXInterface()->write(utils::getAbsoluteFilename(runtimeSettings.workDir, "objective.xml"), objSeriesPI, nEnsemble, JInc3D, "");

	// some final analysis on the objective function terms, contribution by term over all time steps
	multimap<double, string> JMap = multimap<double, string>();
	for (int i=0; i<(int)obj->getNTerm(); i++) {
		double JTerm = 0.0;
		for (int j=0; j<(int)ensembleMap.size(); j++) {
			for (int k=0; k<tsModel->getTimeSeriesTensor()->getNTimeStep(); k++) {
				if (JInc3D[ensembleMap[j]][k][5*i+4]==JInc3D[ensembleMap[j]][k][5*i+4]) {
					JTerm += JInc3D[ensembleMap[j]][k][5*i+4];
				}
			}
		}
		JMap.insert(pair<double,string>(JTerm, obj->getTermIDs()[i]));
	}

	// csv output of objective function summary
	string filename = utils::getAbsoluteFilename(runtimeSettings.workDir, "objective_summary.csv");
	ofstream csvFile(filename.c_str(), ios::out | ios::trunc);

	csvFile << "id,value" << endl;
	for(multimap<double, string>::reverse_iterator p = JMap.rbegin(); p != JMap.rend(); p++) {
		csvFile << p->second << "," << p->first << endl;
	}
	csvFile.close();

	utils::free_dten(JInc3D);
}


void rtcToolsRuntime::reportConstraints()
{
	const string flag_enum = "constraint status = 0,...,15 (combinations by binary AND): 0 = within inactive bounds, 1 = active lower bound, 2 = active upper bound, 4 = violated lower bound, 8 = violated upper bound";

	// optimization variables and bounds
	int n = getN();
	double* x = new double[n];
	double* x_l = new double[n];
	double* x_u = new double[n];

	// constraints and bounds
	int m = getM();
	double* g = new double[m];
	double* g_l = new double[m];
	double* g_u = new double[m];

	// get data
	getInput(n, x);
	simulate(n, x);
	get_bounds_info(n, x_l, x_u, m, g_l, g_u);
	eval_g(n, x, m, g);

	// variables and its bounds
	vector<metadata> var_md(n);
	getN_metadata(var_md);
	for (int i=0; i<n; i++) {
		var_md[i].flag = 0;
		var_md[i].message = "";
		bool violation = false;
		if (x[i] < x_l[i]-runtimeSettings.constraintViolationTol) {
			stringstream ss;
			ss << "lower bound is violated (" << setprecision(12) << x[i] << " < "
			   << setprecision(12) << x_l[i] << " beyond tolerance of " << runtimeSettings.constraintViolationTol << ")";
			var_md[i].message.append(ss.str());
			var_md[i].flag += 4;
			violation = true;
		}
		if (fabs(x[i]-x_l[i]) < runtimeSettings.constraintViolationTol) {
			if (var_md[i].message.size()>0) var_md[i].message.append(", ");
			var_md[i].message.append("lower bound is active constraint");
			var_md[i].flag += 1;
		}
		if (x[i] > x_u[i]+runtimeSettings.constraintViolationTol) {
			stringstream ss;
			ss << "upper bound is violated (" << setprecision(12) << x[i] << " > "
			   << setprecision(12) << x_u[i] << " beyond tolerance of " << runtimeSettings.constraintViolationTol << ")";
			if (var_md[i].message.size()>0) ss << ", ";
			var_md[i].message.append(ss.str());
			var_md[i].flag += 8;
			violation = true;
		}
		if (fabs(x[i]-x_u[i]) < runtimeSettings.constraintViolationTol) {
			if (var_md[i].message.size()>0) var_md[i].message.append(", ");
			var_md[i].message.append("upper bound is active constraint");
			var_md[i].flag += 2;
		}
		if (runtimeSettings.constraintViolationAsError && violation) {
			// report violation to diag.xml
			piDiagInterface::addLine(runtimeSettings.constraintViolationLevel,
				"void rtcToolsRuntime::finish(void) - constraint violation: "
				+ var_md[i].label + ": " + var_md[i].message, RTCTOOLSRUNTIME_CODE);
		}
	}

	// constraints and its bounds
	vector<metadata> con_md(m);
	getM_metadata(con_md);
	for (int i=0; i<m; i++) {
		con_md[i].flag = 0;
		con_md[i].message = "";
		bool violation = false;
		if (g[i] < g_l[i]-runtimeSettings.constraintViolationTol) {
			stringstream ss;
			ss << "lower bound is violated (" << setprecision(12) << g[i] << " < "
			   << setprecision(12) << g_l[i] << " beyond tolerance of " << runtimeSettings.constraintViolationTol << ")";
			con_md[i].message.append(ss.str());
			con_md[i].flag += 4;
			violation = true;
		}
		if (fabs(g[i]-g_l[i]) <= runtimeSettings.constraintViolationTol) {
			if (con_md[i].message.size()>0) con_md[i].message.append(", ");
			con_md[i].message.append("lower bound is active constraint");
			con_md[i].flag += 1;
		}
		if (g[i] > g_u[i]+runtimeSettings.constraintViolationTol) {
			stringstream ss;
			ss << "upper bound is violated (" << setprecision(12) << g[i] << " > "
			   << setprecision(12) << g_u[i] << " beyond tolerance of " << runtimeSettings.constraintViolationTol << ")";
			if (con_md[i].message.size()>0) ss << ", ";
			con_md[i].message.append(ss.str());
			con_md[i].flag += 8;
			violation = true;
		}
		if (fabs(g[i]-g_u[i]) <= runtimeSettings.constraintViolationTol) {
			if (con_md[i].message.size()>0) con_md[i].message.append(", ");
			con_md[i].message.append("upper bound is active constraint");
			con_md[i].flag += 2;
		}
		if (runtimeSettings.constraintViolationAsError && violation) {
			// report violation to diag.xml
			piDiagInterface::addLine(runtimeSettings.constraintViolationLevel,
				"void rtcToolsRuntime::finish(void) - constraint violation: "
				+ con_md[i].label + ": " + con_md[i].message, RTCTOOLSRUNTIME_CODE);
		}
	}

	// count time series for xml output
	int nConSeries = 4*obj->getNVariable() + 5*obj->getNVariableRateOfChange()
		+ 4*obj->getNVariableAverage() + 4*obj->getNState();

	// allocate and initialize memory
	vector<string> conSeriesCSV = vector<string>(nConSeries);
	vector<piTimeSeries> conSeriesPI = vector<piTimeSeries>(nConSeries);
	double*** conData = utils::dten(nEnsemble, tsModel->getTimeSeriesTensor()->getNTimeStep(), nConSeries);
	for (int i=0; i<nEnsemble; i++) {
		for (int j=0; j<tsModel->getTimeSeriesTensor()->getNTimeStep(); j++) {
			for (int k=0; k<nConSeries; k++) {
				conData[i][j][k] = std::numeric_limits<double>::quiet_NaN();
			}
		}
	}
	stringContainer sc;

	// fill data
	int indx = 0;

	for (int i=0; i<obj->getNVariable(); i++) {

		objectiveFunction::variable* v = obj->getVariable(i);

		conSeriesCSV[indx] = v->id + "_x";
		conSeriesCSV[indx+1] = v->id + "_x_L";
		conSeriesCSV[indx+2] = v->id + "_x_U";
		conSeriesCSV[indx+3] = v->id + "_x_flag";

		// meta data of pi time series, if not configures XML element will be empty
		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(v->id);
		if (pits) {

			// if a pi time series definition is available,
			// copy it, add an additional qualifier and set the index
			conSeriesPI[indx] = *pits;
			conSeriesPI[indx].addQualifierID("x");

			conSeriesPI[indx+1] = *pits;
			conSeriesPI[indx+1].addQualifierID("x_L");

			conSeriesPI[indx+2] = *pits;
			conSeriesPI[indx+2].addQualifierID("x_U");

			conSeriesPI[indx+3] = *pits;
			conSeriesPI[indx+3].addQualifierID("x_flag");
			conSeriesPI[indx+3].setComment(flag_enum);
		}

		conSeriesPI[indx].setIndex(indx);
		conSeriesPI[indx+1].setIndex(indx+1);
		conSeriesPI[indx+2].setIndex(indx+2);
		conSeriesPI[indx+3].setIndex(indx+3);

		// fill data
		int nStepSim = iEnd-iStart+1;
		int nStepOpt = nStepSim;
		if (v->vAggregation) {
			nStepOpt = v->vAggregation->getNVariable();
		}

		for (int j=0; j<(int)ensembleMap.size(); j++) {

			if (v->ensembleMode == objectiveFunction::JOINT) {
				// joint (default) use of one trajectory for every branch if nEnsemble>1
				for (int k=0; k<nStepOpt; k++) {

					int t = iStart + k;
					if (v->vAggregation) t = v->vAggregation->getStepSim(k);

					conData[ensembleMap[j]][t][indx] = x[k+v->nStart]/v->scalingFactor;
					if (x_l[k+v->nStart]>-1e+19) conData[ensembleMap[j]][t][indx+1] = x_l[k+v->nStart]/v->scalingFactor;
					if (x_u[k+v->nStart]<+1e+19) conData[ensembleMap[j]][t][indx+2] = x_u[k+v->nStart]/v->scalingFactor;
					conData[ensembleMap[j]][t][indx+3] = var_md[k+v->nStart].flag;
				}

			} else if (v->ensembleMode == objectiveFunction::TREE) {
				// scenario tree mode
				vector<double> branch(nStepSim);
				treeInt->total2tree(j, treeInt->getNStepTotal(), (double*)(x+v->nStart), branch);
				for (int k=0; k<nStepOpt; k++) {
					conData[ensembleMap[j]][iStart+k][indx] = branch[k]/v->scalingFactor;
				}
				treeInt->total2tree(j, treeInt->getNStepTotal(), (double*)(x_l+v->nStart), branch);
				for (int k=0; k<nStepOpt; k++) {
					if (branch[k]>-1e+19) conData[ensembleMap[j]][iStart+k][indx+1] = branch[k]/v->scalingFactor;
				}
				treeInt->total2tree(j, treeInt->getNStepTotal(), (double*)(x_u+v->nStart), branch);
				for (int k=0; k<nStepOpt; k++) {
					if (branch[k]<+1e+19) conData[ensembleMap[j]][iStart+k][indx+2] = branch[k]/v->scalingFactor;
				}
				// TODO consider flag output

			} else if (v->ensembleMode == objectiveFunction::INDEPENDENT) {
				// independent variables on every ensemble branch
				for (int k=0; k<nStepOpt; k++) {

					int t = iStart + k;
					if (v->vAggregation) t = v->vAggregation->getStepSim(k);

					conData[ensembleMap[j]][t][indx] = x[k+j*nStepOpt+v->nStart]/v->scalingFactor;
					if (x_l[k+j*nStepOpt+v->nStart]>-1e+19) conData[ensembleMap[j]][t][indx+1] = x_l[k+j*nStepOpt+v->nStart]/v->scalingFactor;
					if (x_u[k+j*nStepOpt+v->nStart]<+1e+19) conData[ensembleMap[j]][t][indx+2] = x_u[k+j*nStepOpt+v->nStart]/v->scalingFactor;
					conData[ensembleMap[j]][t][indx+3] = var_md[k+j*nStepOpt+v->nStart].flag;
				}
			}
		}
		indx += 4;
	}

	for (int i=0; i<obj->getNVariableRateOfChange(); i++) {

		objectiveFunction::variableRateOfChange* rof = obj->getVariableRateOfChange(i);

		conSeriesCSV[indx] = rof->id + "_dcMin";
		conSeriesCSV[indx+1] = rof->id + "_dcMax";
		conSeriesCSV[indx+2] = rof->id + "_dc_L";
		conSeriesCSV[indx+3] = rof->id + "_dc_U";
		conSeriesCSV[indx+4] = rof->id + "_dc_flag";

		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(rof->variableID);
		if (pits) {
			conSeriesPI[indx] = *pits;
			conSeriesPI[indx].addQualifierID(rof->name);
			conSeriesPI[indx].addQualifierID("dcMin");

			conSeriesPI[indx+1] = *pits;
			conSeriesPI[indx+1].addQualifierID(rof->name);
			conSeriesPI[indx+1].addQualifierID("dcMax");

			conSeriesPI[indx+2] = *pits;
			conSeriesPI[indx+2].addQualifierID(rof->name);
			conSeriesPI[indx+2].addQualifierID("dc_L");

			conSeriesPI[indx+3] = *pits;
			conSeriesPI[indx+3].addQualifierID(rof->name);
			conSeriesPI[indx+3].addQualifierID("dc_U");

			conSeriesPI[indx+4] = *pits;
			conSeriesPI[indx+4].addQualifierID(rof->name);
			conSeriesPI[indx+4].addQualifierID("dc_flag");
			conSeriesPI[indx+4].setComment(flag_enum);
		}

		conSeriesPI[indx].setIndex(indx);
		conSeriesPI[indx+1].setIndex(indx+1);
		conSeriesPI[indx+2].setIndex(indx+2);
		conSeriesPI[indx+3].setIndex(indx+3);
		conSeriesPI[indx+4].setIndex(indx+4);

		// fill data
		for (int j=0; j<(int)rof->rof2DArray.size(); j++) {
			for (int k=0; k<(int)rof->rof2DArray[j].size(); k++) {
				if (rof->rof2DArray[j][k].iStep.size() == 0)
					continue;

				int iStep = rof->rof2DArray[j][k].iStep.back();

				if (conData[j][iStep][indx]!=conData[j][iStep][indx]) {
					// first call
					conData[j][iStep][indx] = g[rof->rof2DArray[j][k].ncIndx];
					conData[j][iStep][indx+1] = g[rof->rof2DArray[j][k].ncIndx];
					conData[j][iStep][indx+4] = con_md[rof->rof2DArray[j][k].ncIndx].flag;
				} else {
					// conData[j][iStep][...] has been already filled before, compute min/max values
					conData[j][iStep][indx] = min(conData[j][iStep][indx], g[rof->rof2DArray[j][k].ncIndx]);
					conData[j][iStep][indx+1] = max(conData[j][iStep][indx+1], g[rof->rof2DArray[j][k].ncIndx]);
					conData[j][iStep][indx+4] = (int)conData[j][iStep][indx+4] & con_md[rof->rof2DArray[j][k].ncIndx].flag;
				}
				if (g_l[rof->rof2DArray[j][k].ncIndx]>-1e+19) conData[j][iStep][indx+2] = g_l[rof->rof2DArray[j][k].ncIndx];
				if (g_u[rof->rof2DArray[j][k].ncIndx]<+1e+19) conData[j][iStep][indx+3] = g_u[rof->rof2DArray[j][k].ncIndx];

				if (con_md[rof->rof2DArray[j][k].ncIndx].message.size()>0) {
					sc.addString(j, iStep, indx, con_md[rof->rof2DArray[j][k].ncIndx].message);
				}
			}
		}
		indx += 5;
	}

	for (int i=0; i<obj->getNVariableAverage(); i++) {

		objectiveFunction::variableAverage* av = obj->getVariableAverage(i);

		conSeriesCSV[indx] = av->id + "_cAv";
		conSeriesCSV[indx+1] = av->id + "_cAv_L";
		conSeriesCSV[indx+2] = av->id + "_cAv_U";
		conSeriesCSV[indx+3] = av->id + "_cAv_flag";

		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(av->variableID);
		if (pits) {
			conSeriesPI[indx] = *pits;
			conSeriesPI[indx].addQualifierID(av->name);
			conSeriesPI[indx].addQualifierID("cAv");

			conSeriesPI[indx+1] = *pits;
			conSeriesPI[indx+1].addQualifierID(av->name);
			conSeriesPI[indx+1].addQualifierID("cAv_L");

			conSeriesPI[indx+2] = *pits;
			conSeriesPI[indx+2].addQualifierID(av->name);
			conSeriesPI[indx+2].addQualifierID("cAv_U");

			conSeriesPI[indx+3] = *pits;
			conSeriesPI[indx+3].addQualifierID(av->name);
			conSeriesPI[indx+3].addQualifierID("cAv_flag");
			conSeriesPI[indx+3].setComment(flag_enum);
		}

		conSeriesPI[indx].setIndex(indx);
		conSeriesPI[indx+1].setIndex(indx+1);
		conSeriesPI[indx+2].setIndex(indx+2);
		conSeriesPI[indx+3].setIndex(indx+3);

		// fill data
		for (int j=0; j<(int)av->av2DArray.size(); j++) {
			for (int k=0; k<(int)av->av2DArray[j].size(); k++) {
				if (av->av2DArray[j][k].iStep.size() == 0)
					continue;

				int iStep = av->av2DArray[j][k].iStep.back();
				conData[j][iStep][indx] = g[av->av2DArray[j][k].ncIndx];
				if (g_l[av->av2DArray[j][k].ncIndx]>-1e+19) conData[j][iStep][indx+1] = g_l[av->av2DArray[j][k].ncIndx];
				if (g_u[av->av2DArray[j][k].ncIndx]<+1e+19) conData[j][iStep][indx+2] = g_u[av->av2DArray[j][k].ncIndx];
				conData[j][iStep][indx] = con_md[av->av2DArray[j][k].ncIndx].flag;

				if (con_md[av->av2DArray[j][k].ncIndx].message.size()>0) {
					sc.setString(j, iStep, indx, con_md[av->av2DArray[j][k].ncIndx].message);
				}
			}
		}
		indx += 4;
	}

	for (int i=0; i<obj->getNState(); i++) {

		objectiveFunction::state* s = obj->getState(i);

		conSeriesCSV[indx] = s->id + "_c";
		conSeriesCSV[indx+1] = s->id + "_c_L";
		conSeriesCSV[indx+2] = s->id + "_c_U";
		conSeriesCSV[indx+3] = s->id + "_c_flag";

		piTimeSeries* pits = tsModel->getPISAXInterface()->getPiTimeSeries(s->stateID);
		if (pits) {
			conSeriesPI[indx] = *pits;
			conSeriesPI[indx].addQualifierID(s->name);
			conSeriesPI[indx].addQualifierID("c");

			conSeriesPI[indx+1] = *pits;
			conSeriesPI[indx+1].addQualifierID(s->name);
			conSeriesPI[indx+1].addQualifierID("c_L");

			conSeriesPI[indx+2] = *pits;
			conSeriesPI[indx+2].addQualifierID(s->name);
			conSeriesPI[indx+2].addQualifierID("c_U");

			conSeriesPI[indx+3] = *pits;
			conSeriesPI[indx+3].addQualifierID(s->name);
			conSeriesPI[indx+3].addQualifierID("c_flag");
			conSeriesPI[indx+3].setComment(flag_enum);
		}

		conSeriesPI[indx].setIndex(indx);
		conSeriesPI[indx+1].setIndex(indx+1);
		conSeriesPI[indx+2].setIndex(indx+2);
		conSeriesPI[indx+3].setIndex(indx+3);

		// fill data
		for (int j=0; j<(int)s->state2DArray.size(); j++) {
			for (int k=0; k<(int)s->state2DArray[j].size(); k++) {

				int iStep = s->state2DArray[j][k].iStepResiduum;
				conData[j][iStep][indx] = g[s->state2DArray[j][k].ncIndx]/s->scalingFactor;
				if (g_l[s->state2DArray[j][k].ncIndx]>-1e+19) conData[j][iStep][indx+1] = g_l[s->state2DArray[j][k].ncIndx]/s->scalingFactor;
				if (g_u[s->state2DArray[j][k].ncIndx]<+1e+19) conData[j][iStep][indx+2] = g_u[s->state2DArray[j][k].ncIndx]/s->scalingFactor;
				conData[j][iStep][indx+3] = con_md[s->state2DArray[j][k].ncIndx].flag;

				if (con_md[s->state2DArray[j][k].ncIndx].message.size()>0) {
					sc.setString(j, iStep, indx, con_md[s->state2DArray[j][k].ncIndx].message);
				}
			}
		}
		indx += 4;
	}

	// csv and pi xml write actions
	if (tsModel->getCsvInterface()) {
		tsModel->getCsvInterface()->writeFiles("constraints", conSeriesCSV, conData);
	}
	tsModel->getPISAXInterface()->write(utils::getAbsoluteFilename(runtimeSettings.workDir, "constraints.xml"), conSeriesPI, nEnsemble, conData, "", sc);

	// free memory
	delete[] x;
	delete[] x_l;
	delete[] x_u;
	delete[] g;
	delete[] g_l;
	delete[] g_u;
}


void rtcToolsRuntime::finish(int timeStep /* =0*/)
{
	char buffer[200];

	clock_t cpu3 = clock();

	try {
		// adjoint output in sensitivity and optimization modes
		bool reportObj = false;
		bool reportCon = false;
		bool reportAdj = false;
		for (int i=0; i<(int)runtimeSettings.modeInfo.size(); i++) {
			if (runtimeSettings.modeInfo[i].mode == FIRSTORDERSENSITIVITY ||
				runtimeSettings.modeInfo[i].mode == OPTIMIZE ||
				runtimeSettings.modeInfo[i].mode == CLOSEDLOOP) {
					reportObj = true;
					reportCon = true;
					reportAdj = true;
			}
			else if (runtimeSettings.modeInfo[i].mode == SIMULATE) {
				if (runtimeSettings.modeInfo[i].executeObjectiveFunction) reportObj = true;
				if (runtimeSettings.modeInfo[i].executeConstraints) reportCon = true;
			}
		}

		if (reportObj) {
			if (runtimeSettings.outputMetaData) reportMetaData();

			// write objectives if any
			if (obj->getNTerm() > 0) {
				reportObjectives();
			}
		}

		if (reportCon) {
			// write constraints if any
			reportConstraints();
		}

		// data validation at model termination
		tsModel->getTimeSeriesTensor()->validate(iForecastTime, false);

		// export time series
		if (timeStep > 0) {
			tsModel->write(reportAdj, timeStep, true);
		} else {
			tsModel->write(reportAdj);
		}
	}
	catch (exception &e)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::finish(void) - error - " + string(e.what()), RTCTOOLSRUNTIME_CODE);
		piDiagInterface::write();
		throw;
	}
	catch (...)
	{
		piDiagInterface::addLine(1, "void rtcToolsRuntime::finish(void) - unknown error in rtcTools", RTCTOOLSRUNTIME_CODE);
		piDiagInterface::write();
		throw runtime_error("void rtcToolsRuntime::finish() - unknown error");
	}

	// some remaining timing statistics
	clock_t cpu4 = clock();

	// cpu1 and cpu2 are evaluated in the main function!
	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (pre-processing) = %d ms", (int)(cpu2-cpu1));
	piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (execution) = %d ms", (int)(cpu3-cpu2));
	piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

	snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (post processing) = %d ms", (int)(cpu4-cpu3));
	piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

	piDiagInterface::write();
	piDiagInterface::freeRtcDiagWriter();
}


void rtcToolsRuntime::getTimeSeriesMatrix(int nTimeStep, int nSeries, double *tsm)
{
	if (!runtimeSettings.limitedMemory) {
        timeseries::timeSeriesMatrix* tsm_temp = dynamic_cast<timeseries::timeSeriesMatrix*>(tsModel->getTimeSeriesTensor()->getTimeSeriesMatrix(0));
		double **tsm_intern = tsm_temp->getValueMatrix();

		for (int i=0; i<nTimeStep; i++) {
			for (int j=0; j<nSeries; j++) {
				tsm[i+nTimeStep*j] = tsm_intern[i][j];
			}
		}
	} else {
		piDiagInterface::addLine(1, "void rtcToolsRuntime::getTimeSeriesMatrix - method not available if limitedMemory option is true", RTCTOOLSRUNTIME_CODE);
	}

}


void rtcToolsRuntime::setTimeSeriesMatrix(int nTimeStep, int nSeries, double *tsm)
{
	if (!runtimeSettings.limitedMemory) {
		timeseries::timeSeriesMatrix* tsMatrix= dynamic_cast<timeSeriesMatrix*>(tsModel->getTimeSeriesTensor()->getTimeSeriesMatrix(0));
		double** tsm_intern = tsMatrix->getValueMatrix();
	    for (int i=0; i<nTimeStep; i++) {
			for (int j=0; j<nSeries; j++) {
				tsm_intern[i][j] = tsm[i+nTimeStep*j];
			}
		}
	} else {
		piDiagInterface::addLine(1, "void rtcToolsRuntime::setTimeSeriesMatrix - method not available if limitedMemory option is true", RTCTOOLSRUNTIME_CODE);
	}
}


bool rtcToolsRuntime::showDebug()
{
	return piDiagInterface::getLogLevel() == 4;
}


void rtcToolsRuntime::setPeriod(simPeriodEnum period)
{
	// set the complete runtime as default
	this->iStart = 1;
	this->iEnd = tsModel->getTimeSeriesTensor()->getNTimeStep()-1;

	// apply a reduction to update or forecast period if applicable
	if (period==UPDATE) {
		this->iEnd = this->iForecastTime;
	} else if (period==FORECAST) {
		this->iStart = this->iForecastTime+1;
	}

	// by default simulate all ensemble members
	ensembleMap = vector<int>(nEnsemble);
	pVec = vector<double>(nEnsemble);
	for (int i=0; i<nEnsemble; i++) {
		ensembleMap[i] = i;
		pVec[i] = 1.0/(double)nEnsemble;
	}
	tsModel->getTimeSeriesTensor()->setEnsembleMap(ensembleMap);

	// initialize optimization problem if available
	if (obj!=0 && this->executeObjectiveFunction) {
		if (treeInt!=0) {
			// initialize the scenario tree if available and mofify the ensemble map if applicable
			treeInt->initialize(tsModel->getTimeSeriesTensor(), this->iStart, this->iEnd);
			ensembleMap = treeInt->getEnsembleMap();
			tsModel->getTimeSeriesTensor()->setEnsembleMap(ensembleMap);
			pVec = treeInt->getProbability();
		}
		obj->initialize(this->iStart, this->iEnd);
	}

	// prepare the simulators for each member
	rtcSim = new rtcToolsSimulator[ensembleMap.size()];
	for (int i=0; i<(int)ensembleMap.size(); i++) {
		rtcSim[i] = rtcToolsSimulator(i, tsModel->getTimeSeriesTensor()->getTimeSeriesMatrix(ensembleMap[i]), schema,
			obj, pVec[i], &runtimeSettings);
	}
}


bool rtcToolsRuntime::execute(int timeStep /* =0 */)
{
	char buffer[200];

	try {
		for (int i=0; i<(int)runtimeSettings.modeInfo.size(); i++) {

			// set flags for the execution of the objective function and constraints
			this->executeObjectiveFunction = runtimeSettings.modeInfo[i].executeObjectiveFunction;
			this->executeConstraints = runtimeSettings.modeInfo[i].executeConstraints;

			// set execution period
			this->setPeriod(runtimeSettings.modeInfo[i].period);

			if (runtimeSettings.modeInfo[i].mode==SIMULATE) {

				clock_t sim_start = clock();
				if (timeStep < 1) {
					J = simulate();         // Perform entire simulation

					snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: objective function value = %lf", J);
					piDiagInterface::addLine(3, string(buffer));

				} else {
					simulate(timeStep); // Simulate until given time step

				}
				clock_t sim = clock()-sim_start;

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (simulation) = %d ms", (int)sim);
				piDiagInterface::addLine(3, string(buffer));

			} else if (runtimeSettings.modeInfo[i].mode==FIRSTORDERSENSITIVITY) {

				clock_t sim_start = clock();
				J = simulate();
				clock_t sim = clock()-sim_start;

				clock_t der_start = clock();
				eval_grad_f(0, 0);
				clock_t der = clock()-der_start;

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: objective function value = %lf", J);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (simulation) = %d ms", (int)sim);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (first-order sensitivity) = %d ms", (int)der);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

			} else if (runtimeSettings.modeInfo[i].mode==FIRSTORDERDERIVATIVECHECK) {

				int n = getN();
				double* x = new double[n];

				getInput(n, x);

				bool have_errors = false;
				if (runtimeSettings.modeInfo[i].derivativeCheck_method==RIDDERS) {
					have_errors = check_grad_f_ridders(n, x, runtimeSettings.modeInfo[i].derivativeCheck_pertubation, runtimeSettings.modeInfo[i].derivativeCheck_tol);
				} else if (runtimeSettings.modeInfo[i].derivativeCheck_method==CENTRAL) {
					have_errors = check_grad_f_central(n, x, runtimeSettings.modeInfo[i].derivativeCheck_pertubation, runtimeSettings.modeInfo[i].derivativeCheck_tol);
				}

				delete[] x;

				if (have_errors && runtimeSettings.modeInfo[i].derivativeCheck_exit_on_error) {
					// Errors occurred in the derivative checker.  Exit now.
					exit(1);
				}

			} else if (runtimeSettings.modeInfo[i].mode==OPTIMIZE ||
					   runtimeSettings.modeInfo[i].mode==POSTPROCESSING) {

				clock_t opt_start = clock();
				try {
					optimizers[i]->optimize();
				} catch (exception& e) {
				piDiagInterface::addLine(1, string("rtcToolsRuntime: error in optimization - ") + e.what(), RTCTOOLSRUNTIME_CODE);
				}
				optimizers[i]->write(runtimeSettings.workDir);
				clock_t opt = clock()-opt_start;

				clock_t sim_start = clock();
				J = simulate();
				clock_t sim = clock()-sim_start;
				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: objective function value = %lf", J);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (model predictive control) = %d ms", (int)opt);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (model predictive control, embedded simulation) = %d ms", (int)sim);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);

			} else if (runtimeSettings.modeInfo[i].mode==CLOSEDLOOP) {

				clock_t opt_start = clock();
				int nStep = tsModel->getTimeSeriesTensor()->getNTimeStep()-1;
				int fh = runtimeSettings.forecastHorizon;
				int rh = runtimeSettings.recedingHorizon;
				for (int j=1; j<=nStep-fh+1; j+=rh) {

					this->iStart = j;
					this->iEnd = j+fh-1;

					// optimization
					if (treeInt!=0) treeInt->initialize(tsModel->getTimeSeriesTensor(), this->iStart, this->iEnd);
					obj->initialize(iStart, iEnd);
					optimizers[i]->optimize();

					// simulation
					for (int k=this->iStart; k<this->iStart+rh-1; k++) {
						simulate(k);
					}
				}
				clock_t opt = clock()-opt_start;

				snprintf(buffer, sizeof(buffer), "rtcToolsRuntime: CPU time (closed-loop simulation) = %d ms", (int)opt);
				piDiagInterface::addLine(3, string(buffer), RTCTOOLSRUNTIME_CODE);
			}
		}
	} catch (exception &e) {
		piDiagInterface::addLine(1, "rtcToolsRuntime::execute() - error - " + string(e.what()) + ".", RTCTOOLSRUNTIME_CODE);
		throw;
	} catch (...) {
		piDiagInterface::addLine(1, "rtcToolsRuntime::execute() - unknown error.", RTCTOOLSRUNTIME_CODE);
		throw;
	}


	// returns false for correct execution, true if an error occured
    return piDiagInterface::getIncludeError();
}


void rtcToolsRuntime::assertFile(string filename)
{
	if (utils::fileAvailable(filename)==false) {
		string m = "file not found: " + filename;
		piDiagInterface::addLine(1, "void rtcToolsRuntime::assertFile(string filename), " + filename, RTCTOOLSRUNTIME_CODE);
		throw runtime_error(m.c_str());
	}
}


long long rtcToolsRuntime::getTime(fews::DateTimeComplexType::DateType date, fews::DateTimeComplexType::TimeType time)
{
	int t[7];

	t[0] = date.year();
	t[1] = date.month();
	t[2] = date.day();
	t[3] = time.hours();
	t[4] = time.minutes();
	t[5] = (int)time.seconds();
	t[6] = 0;

	return utils::date2time(t);
}


long long rtcToolsRuntime::getTimeStep(fews::TimeStepComplexType::UnitType unit)
{
	long long dt = -1;

	if (unit.compare("week")==0) {
		dt = 7*24*3600*1000;
	} else if (unit.compare("day")==0) {
		dt = 24*3600*1000;
	} else if (unit.compare("hour")==0) {
		dt = 3600*1000;
	} else if (unit.compare("minute")==0) {
		dt = 60*1000;
	} else if (unit.compare("second")==0) {
		dt = 1000;
	} else {
		throw runtime_error("error - time step unit not implemented in long long rtcToolsHelperFunctions::getTimeStep(TimeStepComplexType::UnitType unit)");
	}

	return dt;
}