/*!
\file
\brief Robust
\author Vasek Smidl 

 */

#include "estim/arx.h"
#include "robustlib.h"
#include <vector>
#include <iostream>
#include <fstream>
//#include <itpp/itsignal.h>
#include "windows.h"
#include "ddeml.h"
#include "stdio.h"
#include <dirent.h>
//#include <itpp/itoptim.h>

//#include "DDEClient.h"
//#include <conio.h>


using namespace itpp;
using namespace bdm;

//const int emlig_size = 2;
//const int utility_constant = 5;

const int max_model_order         = 2;
const double apriorno             = 0.01;
const int max_window_size         = 40;
const int utility_order           = 25;
const int prediction_time         = 30;
const double min_utility_argument = 0.001;
const double max_investment       = 3.0;
const int sample_size			  = 5000;
const char* commodity		      = "TY\\";

/*
HDDEDATA CALLBACK DdeCallback(
	UINT uType,     // Transaction type.
	UINT uFmt,      // Clipboard data format.
	HCONV hconv,    // Handle to the conversation.
	HSZ hsz1,       // Handle to a string.
	HSZ hsz2,       // Handle to a string.
	HDDEDATA hdata, // Handle to a global memory object.
	DWORD dwData1,  // Transaction-specific data.
	DWORD dwData2)  // Transaction-specific data.
{
	return 0;
}

void DDERequest(DWORD idInst, HCONV hConv, char* szItem)
{
	HSZ hszItem = DdeCreateStringHandle(idInst, szItem, 0);
	HDDEDATA hData = DdeClientTransaction(NULL,0,hConv,hszItem,CF_TEXT, 
									XTYP_ADVSTART,TIMEOUT_ASYNC , NULL); //TIMEOUT_ASYNC
	if (hData==NULL)
	{
		printf("Request failed: %s\n", szItem);
	}

	if (hData==0)
	{
		printf("Request failed: %s\n", szItem);
	}
}

DWORD WINAPI ThrdFunc( LPVOID n )
{    
    return 0;
}
*/
vector<char*> listFiles(char* dir){
        vector<char*> files;
		DIR *pDIR;
        struct dirent *entry;
        if( pDIR=opendir(dir) ){
                while(entry = readdir(pDIR)){
                        if( strcmp(entry->d_name, ".") != 0 && strcmp(entry->d_name, "..") != 0 )
						{							
							char *file_string = new char[strlen(entry->d_name) + 1];
							strcpy(file_string, entry->d_name);												
							files.push_back(file_string);
						}
                }
                closedir(pDIR);
        }
		return files;
}

double valueCRRAUtility(const double &position, const pair<vec,vec> &samples, const int order)
{
	double value = 0;
	set<double> values;
			
	for(int i=0;i<samples.first.length();i++)
	{
		double sample      = samples.second.get(i);
		double probability = samples.first.get(i);
		if((position*sample+1)>min_utility_argument)
		{
			values.insert(probability*sample/pow(position*sample+1,order+1));
		}
		else
		{
			values.insert(probability*(min_utility_argument-1)/position/pow(min_utility_argument,order+1));
		}
	}	
	
	for(set<double>::iterator val_ref = values.begin();val_ref!=values.end();val_ref++)
	{
		value+=(*val_ref);
	}

	return value;
}

double gradientCRRAUtility(const double &position, const pair<vec,vec> &samples, const int order)
{	
	double value = 0;
	set<double> values;

	for(int i=0;i<samples.first.length();i++)
	{
		double sample      = samples.second.get(i);
		double probability = samples.first.get(i);
		if((position*sample+1)>min_utility_argument)
		{
			values.insert((-(order+1)*pow(sample,2)*probability)/pow(position*sample+1,order+2));
		}
	}

	for(set<double>::iterator val_ref = values.begin();val_ref!=values.end();val_ref++)
	{
		value+=(*val_ref);
	}

	return value;
}

double newtonRaphson(double startingPoint, double epsilon, pair<vec,vec> samples, int order)
{
	/*
	if(samples.length()>800)
	{
		samples.del(801,samples.size()-1);
	}
	*/
	
	int count = 0;
	
	bool epsilon_reached = false;

	while(count<1000&&!epsilon_reached)
	{
		double cur_value = valueCRRAUtility(startingPoint,samples,order);
		double cur_gradient = gradientCRRAUtility(startingPoint,samples,order);

		startingPoint = startingPoint - cur_value/cur_gradient;

		if(cur_value<epsilon)
		{
			epsilon_reached = true;
		}
	}

	if(count==100)
	{
		return startingPoint; // can be different!
	}
	else
	{
		return startingPoint;
	}
}

class model
{
public:
	set<pair<int,int>> ar_components;

	// Best thing would be to inherit the two models from a single souce, this is planned, but now structurally
	// problematic.
	RARX* my_rarx; //vzmenovane parametre pre triedu model
	ARXwin* my_arx;

	bool has_constant;
	int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit
	int  predicted_channel;
	mat* data_matrix;
	vec  predictions;
	char name[80];

	double previous_lognc;
	
	model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor
		  bool robust, 
		  bool has_constant, 
		  int window_size, 
		  int predicted_channel,
		  mat* data_matrix)
	{		
		this->ar_components.insert(ar_components.begin(),ar_components.end());

		strcpy(name,"M");

		for(set<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++)
		{
			char buffer1[2];
			char buffer2[2];
			itoa((*ar_ref).first,buffer1,10);
			itoa((*ar_ref).second,buffer2,10);

			strcat(name,buffer1);
			strcat(name,buffer2);
			strcat(name,"_");
		}

		this->has_constant      = has_constant;

		if(has_constant)
		{
			strcat(name,"C");
		}

		this->window_size       = window_size;
		this->predicted_channel = predicted_channel;
		this->data_matrix       = data_matrix;

		if(robust)
		{
			previous_lognc = 0;
			strcat(name,"R");

			if(has_constant)
			{
				my_rarx = new RARX(ar_components.size()+1,window_size,true,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+3);
				my_arx  = NULL;
			}
		else
			{
				my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+2);
				my_arx  = NULL;
			}
		}
		else
		{
			my_rarx = NULL;
			my_arx  = new ARXwin();
			mat V0;

			if(has_constant)
			{				
				V0  = apriorno * eye(ar_components.size()+2); //aj tu konst
				V0(0,0) = 0;
				my_arx->set_constant(true);
				my_arx->set_statistics(1, V0, V0.rows()+1); 
			}
			else
			{
				
				V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
				V0(0,0) = 0;
				//V0(0,1) = -0.01;
				//V0(1,0) = -0.01;
				my_arx->set_constant(false);
				my_arx->set_statistics(1, V0, V0.rows()+1); 
				
			}						
			
			my_arx->set_parameters(window_size);
			my_arx->validate();

			previous_lognc = my_arx->posterior().lognc();

			/*
			vec mean = my_arx->posterior().mean();
			cout << mean << endl;
			*/
		}
	}

	~model()
	{
		if(my_rarx!=NULL)
		{
			delete my_rarx;
		}
		else
		{
			delete my_arx;
		}
	}

	void data_update(int time) 
	{
		vec data_vector;
		for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
		{  
			data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
		}
		if(my_rarx!=NULL)
		{	
			data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
			my_rarx->bayes(data_vector);
		}
		else
		{
			vec pred_vec;
			pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
			my_arx->bayes(pred_vec,data_vector);
		}
	}

	pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to
	{
		vec condition_vector;
		for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
		{
			condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1));
		}

		if(my_rarx!=NULL)
		{
			pair<vec,mat> imp_samples = my_rarx->posterior->sample(sample_size,false);

			//cout << imp_samples.first << endl;			
			
			vec sample_prediction;			
			for(int t = 0;t<imp_samples.second.cols();t++)
			{
				vec lap_sample = condition_vector;
				
				if(has_constant)
				{
					lap_sample.ins(lap_sample.size(),1.0);
				}
				
				lap_sample.ins(lap_sample.size(),(*LapRNG)());

				sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));				
			}

			return pair<vec,vec>(imp_samples.first,sample_prediction);
		}	
		else
		{
			mat samples = my_arx->posterior().sample_mat(sample_size);			
			
			vec sample_prediction;
			for(int t = 0;t<sample_size;t++)
			{
				vec gau_sample = condition_vector;
				
				if(has_constant)
				{
					gau_sample.ins(gau_sample.size(),1.0);
				}
				
				gau_sample.ins(gau_sample.size(),randn());

				vec param_sample = samples.get_col(t);
				param_sample.set(param_sample.size()-1,sqrt(param_sample[param_sample.size()-1]));
				sample_prediction.ins(0,gau_sample*param_sample);				
			}

			return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction);
		}
	
	}


	static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels)
	{
		set<set<pair<int,int>>> created_model_types;		

		if(max_order == 1)//ukoncovacia vetva
		{			
			for(int channel = 0;channel<number_of_channels;channel++)//pre AR 1 model vytvori kombinace kanalov v prvom kroku poyadu
			{
				set<pair<int,int>> returned_type;
				returned_type.insert(pair<int,int>(channel,1));  //??
				created_model_types.insert(returned_type);
			}

			return created_model_types;
		}
		else
		{
			created_model_types = possible_models_recurse(max_order-1,number_of_channels);//tu uz mame ulozene kombinace o jeden krok dozadu  //rekuryivne volanie
			set<set<pair<int,int>>> returned_types;
			
			for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++)
			{				
				
				for(int order = 1; order<=max_order; order++)
				{
					for(int channel = 0;channel<number_of_channels;channel++)
					{
						set<pair<int,int>> returned_type;
						pair<int,int> new_pair = pair<int,int>(channel,order);//??
						if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) //??
						{
							returned_type.insert((*model_ref).begin(),(*model_ref).end()); //co vlozi na zaciatok retuned_type?
							returned_type.insert(new_pair);							

							returned_types.insert(returned_type);							
						}
					}
				}
			}

			created_model_types.insert(returned_types.begin(),returned_types.end());

			return created_model_types;
		}		
	}
};

// ****************************************************
//			     MAIN MAIN MAIN MAIN MAIN
// ****************************************************
int main ( int argc, char* argv[] ) 
{
	
	// EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from 
    // y_t=0.95*y_(t-1)+0.05*y_(t-2)+0.2*e_t, where e_t is normally, student(4) and cauchy distributed. It
    // can be compared to the classical setup.
	
	itpp::Laplace_RNG LapRNG = Laplace_RNG();

	char* folder_string  =  "C:\\RobustExperiments\\"; // "C:\\dataADClosePercDiff"; //
	char* data_folder    = "data\\";
	char* results_folder = "results\\";	

	char dfstring[150];
	strcpy(dfstring,folder_string);
	strcat(dfstring,data_folder);
	strcat(dfstring,commodity);
	vector<char*> files = listFiles(dfstring);
	
	for(int contract=0;contract<files.size();contract++)
	{
		char *cdf_str = new char[strlen(dfstring) + 1];
		strcpy(cdf_str, dfstring);		
		strcat(cdf_str,files[contract]);	
	
		mat data_matrix;
		ifstream myfile(cdf_str);
		if (myfile.is_open())
		{		
			string line;
			while(getline(myfile,line))
			{			
				vec data_vector;
				while(line.find(',') != string::npos) //zmenil som ciarku za medzeru
				{
					//line.erase(0,1); // toto som sem pridal
					int loc2 = line.find('\n');
					int loc  = line.find(',');
					data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));				
					line.erase(0,loc+1);					
				}

				data_matrix.ins_row(data_matrix.rows(),data_vector);
			}		

			myfile.close();	
		}
		else
		{
			cout << "Can't open data file!" << endl;
		}
	
		//konec nacitavania dat
		set<set<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,data_matrix.rows()); //volanie funkce kde robi kombinace modelov
													//to priradime do model_types, data_matrix.row urcuje pocet kanalov dat
		vector<model*> models;
		for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++)
		{// prechadza rozne typy kanalov, a poctu regresorov
			for(int window_size = max_window_size-1;window_size < max_window_size;window_size++)
			{
				//if(model_type->size()<max_model_order)
				//{
					models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy
					models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));
				//}
				models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix));
				models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));		
			}

			//set<pair<int,int>> empty_list;
			//models.push_back(new model(empty_list,false,true,100,0,&data_matrix));
		}
	
		mat result_lognc;
		// mat result_preds;

		ofstream myfilew;
		char rfstring[150];					
		strcpy(rfstring,folder_string);
		strcat(rfstring,results_folder);
		strcat(rfstring,commodity);
		strcat(rfstring,files[contract]);

		/*
		char predstring[150];					
		strcpy(predstring,folder_string);
		strcat(predstring,results_folder);
		strcat(predstring,commodity);
		strcat(predstring,"PRED");
		strcat(predstring,files[contract]);
		*/

		for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols() 
		{ 	
			cout << "Steps: " << time << endl;
		
			if(time == max_model_order)
			{ 
				/*
				myfilew.open(rfstring,ios::app);
				for(int i = 0;models.size();i++)
				{
					for(set<pair<int,int>>::iterator ar_ref = models[i]->ar_components.begin();ar_ref != models[i]->ar_components.end();ar_ref++)
					{
						myfilew << (*ar_ref).second << (*ar_ref).first;							
					}

					myfilew << ".";

					if(models[i]->my_arx == NULL)
					{
						myfilew << "31";
					}
					else
					{
						myfilew << "30";
					}
					
					if(models[i]->has_constant)
					{
						myfilew << "61";
					}
					else
					{
						myfilew << "60";
					}

					myfilew << ",999,999,999,";
				}

				myfilew << "888" << endl;
				myfilew.close();
				*/
			}
		
			//pocet stlpcov data_matrix je pocet casovych krokov
			double cur_loglikelihood;
			// vec preds;
			int  prev_samples_nr;
			bool previous_switch = true;
			for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++) //.begin()+1;model_ref++)
			{//posuvam s apo models, co je pole modelov urobene o cyklus vyssie. Teda som v case time a robim to tam pre vsetky typy modelov, kombinace regresorov
				(*model_ref)->data_update(time); //pozret sa preco je toto tu nutne

				//cout << "Updated." << endl;				

				if((*model_ref)->my_rarx!=NULL)	//vklada normalizacni faktor do cur_res_lognc
				{
					//cout << "Maxlik vertex(robust):" << (*model_ref)->my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
					cur_loglikelihood = (*model_ref)->my_rarx->posterior->_ll();				
				}
				else
				{
					//cout << "Mean(classical):" << (*model_ref)->my_arx->posterior().mean() << endl;
					double cur_lognc = (*model_ref)->my_arx->posterior().lognc();
					cur_loglikelihood = cur_lognc-(*model_ref)->previous_lognc;
								
					(*model_ref)->previous_lognc = cur_lognc;				
				}
				
				/*
			if(time == max_window_size-1)
			{
				//***********************
				int sample_size = 100000;
				//***********************
				
				pair<vec,mat> samples;
				if((*model_ref)->my_arx!=NULL)
				{
					mat samp_mat = (*model_ref)->my_arx->posterior().sample_mat(sample_size);
					samples = pair<vec,mat>(ones(samp_mat.cols()),samp_mat);
				}
				else
				{
					samples = (*model_ref)->my_rarx->posterior->sample(sample_size,true);					
				}

				char fstring[80];					
				strcpy(fstring,file_string);
				strcat(fstring,(*model_ref)->name);
				strcat(fstring,".txt");
				
				//cout << samples.first << endl;
				
				myfilew.open(fstring,ios::app);
				
				
				//for(int i = 0;i<samples.first.size();i++)
				//{
				//	myfilew << samples.first.get(i) << ",";
				//}
				//myfilew << endl;
				

				for(int j = 0;j<samples.second.rows()+1;j++)
				{			
					for(int i = 0;i<samples.second.cols();i++)
					{
						if(j!=samples.second.rows())
						{
							myfilew << samples.second.get(j,i) << ",";
						}
						
						//else
						//{
						//	myfilew << "0,";
						//}
						
					}				
					myfilew << endl;
				}

				cout << "*************************************" << endl;

				myfilew.close();				
			}
			*/				

				if(time>prediction_time)
				{			
					int samples_nr;					
					if(previous_switch)
					{
						samples_nr = sample_size;
					}
					else
					{
						samples_nr = prev_samples_nr;
					}
					
					
					// PREDICTIONS
					pair<vec,vec> predictions = (*model_ref)->predict(samples_nr,time,&LapRNG);

					/*
					myfilew.open(predstring,ios::app);
					for(int i=0;i<10000;i++)
					{
						if(i<predictions.second.size())
						{
							myfilew << predictions.second.get(i) << ",";
						}
						else
						{
							myfilew << ",";
						}
					}
					myfilew << endl;
					myfilew.close();
					*/

					if(previous_switch)
					{
						prev_samples_nr = predictions.second.size();
						samples_nr = prev_samples_nr;
					}

					previous_switch = !previous_switch;

					double optimalInvestment = newtonRaphson(0,0.00001,predictions,utility_order);

					if(abs(optimalInvestment)>max_investment)
					{
						optimalInvestment = max_investment*sign(optimalInvestment);
					}


					/*
					vec utilityValues;
					for(int j=0;j<1000;j++)
					{
						utilityValues.ins(utilityValues.length(),valueCRRAUtility(-0.5+0.001*j, predictions.second, utility_order));
					}*/

					double avg_prediction = (predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()));

					(*model_ref)->predictions.ins((*model_ref)->predictions.size(),avg_prediction);
						
					myfilew.open(rfstring,ios::app);
				
					/*
					for(int j=0;j<utilityValues.length();j++)
					{
						myfilew << utilityValues.get(j) << ",";
					}
					myfilew << endl;
					*/

					myfilew << avg_prediction << "," << optimalInvestment << "," << samples_nr << "," << cur_loglikelihood << ","; 
					myfilew.close();
				}			
			}

			if(time>prediction_time&&(time+1)<data_matrix.cols())
			{
				// REAL PRICE		
				myfilew.open(rfstring,ios::app);
				myfilew << data_matrix.get(0,time+1) << endl;
				myfilew.close();
			}	
		}

		for(vector<model*>::reverse_iterator model_ref = models.rbegin();model_ref!=models.rend();model_ref++)
		{
			delete *model_ref;
		}
	}

	return 0;
}