/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//  .h file          ////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


#pragma once


#define MKL_DIR "C:/Progra~1/Intel/COMPOS~3/MKL/"

#include "C:/Progra~1/Intel/COMPOS~3/MKL/include/mkl_types.h"
#include "C:/Progra~1/Intel/COMPOS~3/MKL/include/mkl_rci.h"


#ifdef _WIN64
	#pragma comment(lib, MKL_DIR "lib/intel64/mkl_intel_c_dll")
	#pragma comment(lib, MKL_DIR "lib/intel64/mkl_intel_thread_dll")
	#pragma comment(lib, MKL_DIR "lib/intel64/mkl_core_dll")
	//#pragma comment(lib, MKL_DIR "lib/intel64/libiomp5md")
#else
	#pragma comment(lib, MKL_DIR "lib/ia32/mkl_intel_c_dll")
	#pragma comment(lib, MKL_DIR "lib/ia32/mkl_intel_thread_dll")
	#pragma comment(lib, MKL_DIR "lib/ia32/mkl_core_dll")
//	#pragma comment(lib, "C:/Progra~1/Intel/COMPOS~3/redist/ia32/compiler/libiomp5md")
#endif

//var for random number generator
static unsigned int nextSeed = 1;

class Test
{
public:
	Test(void);
	~Test(void);
};



/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//  .cpp file        ////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include "Test.h"

#include <iostream>
#include <fstream>
#include <cstring>
#include <cmath>
#include <cassert>
#include <cstdlib>


//Custom random number generator
//from POSIX.1-2001
void CustomRandSet(unsigned int seed)
{
	nextSeed = seed;
}

int CustomRand()
{
	nextSeed = nextSeed*1103515245 + 12345;
	return int((nextSeed/65536)%RAND_MAX);
}

double randomDouble(double lowerBound, double upperBound)
{
	return (CustomRand() * (upperBound - lowerBound)) / RAND_MAX + lowerBound;
}

//configuration parameters to the non-linear solvers
struct MKLNonLinearConfig
{
	double epsilon[6];
	MKL_INT maxIterations;
	MKL_INT maxTrialIterations;

	double initialStepBound;
	double jacobianPrecision;
};

//output data from the non-linear solvers
struct MKLNonLinearResults
{
	MKL_INT iterationsUsed;
	MKL_INT stopCriterion;
	double initialResiduals;
	double finalResiduals;
};

//data used inside the objective function
struct NonLinearData
{
	double targetAValue;
	double targetBValue;
	double lowerXBounds;
	double upperXBounds;
};

//a reasonably simple non-linear function:  a and b are parameters,  solving
//for x, the position at which to evaluate the funtion
double NonLinearFunction(double a, double b, double x)
{
	return a * std::exp(-b * x);
}

//This objective function samples a number of points between the bounds 
//specified in the user data, and gives the residuals for the solver's current 
//param compared with the target param
void MKLObjectiveFunction(MKL_INT* pNumResiduals, MKL_INT* pNumParams, double* params, double* residuals, void* userData)
{
	const NonLinearData* data = static_cast<NonLinearData*>(userData);
	
	const double stepIncrement = (data->upperXBounds - data->lowerXBounds) / ((*pNumResiduals)-1);

	assert(*pNumParams == 2);

	double currentX = data->lowerXBounds;

	for (size_t i = 0; i < size_t(*pNumResiduals); ++i)
	{
		residuals[i] = NonLinearFunction(data->targetAValue, data->targetBValue, currentX) 
			- NonLinearFunction(params[0], params[1], currentX);
		currentX += stepIncrement;
	}
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//solver without bounds/constraints
void MKLNonLinearSolve(MKL_INT numParams, MKL_INT numResiduals, MKLNonLinearConfig* config, NonLinearData* userData, double* inoutParams, MKLNonLinearResults* results)
{
	double* residuals = new double[numResiduals];
	double* jacobian = new double[numParams*numResiduals];

	_TRNSP_HANDLE_t mklOperationHandle = 0;

	MKL_INT initResult = dtrnlsp_init(&mklOperationHandle, &numParams, &numResiduals, inoutParams, config->epsilon, &config->maxIterations, &config->maxTrialIterations, &config->initialStepBound);
	assert(initResult == TR_SUCCESS);

	MKL_INT rciRequest = 0;

	do
	{
		MKL_INT solveResult = dtrnlsp_solve(&mklOperationHandle, residuals, jacobian, &rciRequest);
		assert(solveResult == TR_SUCCESS);

		switch(rciRequest)
		{
		case 1:	
			MKLObjectiveFunction(&numResiduals, &numParams, inoutParams, residuals, static_cast<void*>(userData));
			break;
		case 2:
			MKL_INT jacobianResult = djacobix(&MKLObjectiveFunction, &numParams, &numResiduals, jacobian, inoutParams, &config->jacobianPrecision, static_cast<void*>(userData));
			assert(jacobianResult == TR_SUCCESS);
			break;
		}
	} while(rciRequest >= 0);

		

	MKL_INT getResult = dtrnlsp_get(&mklOperationHandle, &results->iterationsUsed, &results->stopCriterion, &results->initialResiduals, &results->finalResiduals);
	assert(getResult == TR_SUCCESS);

	MKL_INT deleteResult = dtrnlsp_delete(&mklOperationHandle);
	assert(deleteResult = TR_SUCCESS);
	
	delete[] residuals;
	residuals = 0;
	delete[] jacobian;
	jacobian = 0;
}



///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//solver without bounds/constraints
void MKLNonLinearSolveWithConstraints(MKL_INT numParams, MKL_INT numResiduals, MKLNonLinearConfig* config, NonLinearData* userData, double* lowerBounds, double* upperBounds, double* inoutParams, MKLNonLinearResults* results)
{
	double* residuals = new double[numResiduals];
	double* jacobian = new double[numParams*numResiduals];

	_TRNSP_HANDLE_t mklOperationHandle = 0;

	MKL_INT initResult = dtrnlspbc_init(&mklOperationHandle, &numParams, &numResiduals, inoutParams, lowerBounds, upperBounds, config->epsilon, &config->maxIterations, &config->maxTrialIterations, &config->initialStepBound);
	assert(initResult == TR_SUCCESS);

	
	MKL_INT rciRequest = 0;

	do
	{
		MKL_INT solveResult = dtrnlspbc_solve(&mklOperationHandle, residuals, jacobian, &rciRequest);
		assert(solveResult == TR_SUCCESS);

		switch(rciRequest)
		{
		case 1:	
			MKLObjectiveFunction(&numResiduals, &numParams, inoutParams, residuals, static_cast<void*>(userData));
			break;
		case 2:
			MKL_INT jacobianResult = djacobix(&MKLObjectiveFunction, &numParams, &numResiduals, jacobian, inoutParams, &config->jacobianPrecision, static_cast<void*>(userData));
			assert(jacobianResult == TR_SUCCESS);
			break;
		}
	} while(rciRequest >= 0);


	MKL_INT getResult = dtrnlspbc_get(&mklOperationHandle, &results->iterationsUsed, &results->stopCriterion, &results->initialResiduals, &results->finalResiduals);
	assert(getResult == TR_SUCCESS);

	MKL_INT deleteResult = dtrnlspbc_delete(&mklOperationHandle);
	assert(deleteResult = TR_SUCCESS);
	
	delete[] residuals;
	residuals = 0;
	delete[] jacobian;
	jacobian = 0;
}



/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Called by main

Test::Test(void)
{
	//set up the unbounded and bounded tests
 
//	unsigned int randomSeed = 1;			//   worked
//	unsigned int randomSeed = 2;			//   worked 
//	unsigned int randomSeed = 3;			// error
//	unsigned int randomSeed = 13;			//   worked
//	unsigned int randomSeed = 20;			// error
	unsigned int randomSeed = 1271347576;	// error
//	unsigned int randomSeed = 2945621217;	// error
//	unsigned int randomSeed = 1945621217;	// error
//	unsigned int randomSeed = 11111111111;	//   worked 
//	unsigned int randomSeed = 29458621217;  // error

	
	CustomRandSet(randomSeed);
	std::cout << "Seed = " << randomSeed << std::endl;

	MKLNonLinearConfig config;
	config.maxIterations = 1000;
	config.maxTrialIterations  = 1000;
	for (int i = 0; i < 6; ++i)
	{
		config.epsilon[i] = 0.00001;
	}
	config.initialStepBound = 100.0;
	config.jacobianPrecision = 0.00001; 

	//specify data to be used in objective function
	NonLinearData userData;
	userData.targetAValue = randomDouble(-10.0, 10.0);
	userData.targetBValue = randomDouble(-3.0, 3.0);
	userData.lowerXBounds = -3.0;
	userData.upperXBounds = 3.0;

	//set up starting parameters based on target data - note a copy is kept for comparison at the end
	double startingParameters[2];
	double unconstrainedSolverResults[2];
	double constrainedSolverResults[2];
		
	startingParameters[0] = userData.targetAValue*randomDouble(0.5, 1.5);
	startingParameters[1] = userData.targetBValue*randomDouble(0.5, 1.5);

	//place loose constaints on the parameters
	double lowerBounds[2] = {-100.0, -100.0};
	double upperBounds[2] = {100.0, 100.0};

	MKLNonLinearResults unconstrainedResults;
	MKLNonLinearResults constrainedResults;

	//solve with constraints
	memcpy(constrainedSolverResults, startingParameters, sizeof(double)*2);
	MKLNonLinearSolveWithConstraints(2, 10, &config, &userData, lowerBounds, upperBounds, constrainedSolverResults, &constrainedResults);

	//solve with no constraints 
	memcpy(unconstrainedSolverResults, startingParameters, sizeof(double)*2);
	MKLNonLinearSolve(2, 10, &config, &userData, unconstrainedSolverResults, &unconstrainedResults);

	//print out results 
	std::cout << std::endl << std::endl;
	std::cout << "--------------- Unconstrained Solver ---------------" << std::endl;
	std::cout << "Initial Parameters: a = " << startingParameters[0] << "\tb = " << startingParameters[1] << std::endl;
	std::cout << "Target Parameters:  a = " << userData.targetAValue << "\tb = " << userData.targetBValue << std::endl;
	std::cout << "Final Parameters:   a = " << unconstrainedSolverResults[0] << "\tb = " << unconstrainedSolverResults[1] << std::endl;
//	std::cout << "Results: " << std::endl;
//	std::cout << "\tIterations =        " << unconstrainedResults.iterationsUsed << std::endl;
	std::cout << "\tStop Condition =    " << unconstrainedResults.stopCriterion << std::endl;
//	std::cout << "\tInitial Residuals = " << unconstrainedResults.initialResiduals << std::endl;
	std::cout << "\tFinal Residuals =   " << unconstrainedResults.finalResiduals << std::endl;
	std::cout << std::endl << std::endl;

	std::cout << "--------------- Constrained Solver -----------------" << std::endl;
	std::cout << "Initial Parameters: a = " << startingParameters[0] << "\tb = " << startingParameters[1] << std::endl;
	std::cout << "Target Parameters:  a = " << userData.targetAValue << "\tb = " << userData.targetBValue << std::endl;
	std::cout << "Final Parameters:   a = " << constrainedSolverResults[0] << "\tb = " << constrainedSolverResults[1] << std::endl;
//	std::cout << "Results: " << std::endl;
//	std::cout << "\tIterations =        " << constrainedResults.iterationsUsed << std::endl;
	std::cout << "\tStop Condition =    " << constrainedResults.stopCriterion << std::endl;
//	std::cout << "\tInitial Residuals = " << constrainedResults.initialResiduals << std::endl;
	std::cout << "\tFinal Residuals =   " << constrainedResults.finalResiduals << std::endl;
	std::cout << std::endl << std::endl;
}

Test::~Test(void)
{
}





int main(int argc, char* argv[])
{
	Test test;

	return 0;
}