/*!
\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 "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.001;

/*
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;
}
*/

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;
	
	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());
		this->has_constant      = has_constant;
		this->window_size       = window_size;
		this->predicted_channel = predicted_channel;
		this->data_matrix       = data_matrix;

		if(robust)
		{
			if(has_constant)
			{
				my_rarx = new RARX(ar_components.size()+1,window_size,true,apriorno*sqrt(2.0),apriorno*sqrt(2.0),ar_components.size()+4);
				my_arx  = NULL;
			}
		else
			{
				my_rarx = new RARX(ar_components.size(),window_size,false,apriorno*sqrt(2.0),apriorno*sqrt(2.0),ar_components.size()+3);
				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);				
			}
			else
			{
				
				V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
				//V0(0,0) = 0;
				my_arx->set_constant(false);				
				
			}

			my_arx->set_statistics(1, V0, V0.rows()+2); 			
			my_arx->set_parameters(window_size);
			my_arx->validate();
		}
	}

	void data_update(int time) //vlozime cas a ono vlozi do data_vector podmineky(conditions) a predikce, ktore pouzije do bayes
	{
		vec data_vector;
		for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
		{  //ar?iterator ide len od 1 pod 2, alebo niekedy len 1
			data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
// do data vector vlozi pre dany typ regresoru prislusne cisla z data_matrix. Ale ako? preco time-ar_iterator->second?
		}
		if(my_rarx!=NULL)
		{	//pre robust priradi az tu do data_vector aj predikciu
			data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
			my_rarx->bayes(data_vector);
		}
		else
		{
			vec pred_vec;//tu sa predikcia zadava zvlast
			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->importance_sample(sample_size);

			//cout << imp_samples.first << endl;
			
			vec sample_prediction;			
			for(int t = 0;t<sample_size;t++)
			{
				vec lap_sample = condition_vector;
				
				if(has_constant)
				{
					lap_sample.ins(lap_sample.size(),1.0);
				}
				
				lap_sample.ins(0,(*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());

				sample_prediction.ins(0,gau_sample*samples.get_col(t));				
			}

			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;
		}		
	}
};




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();	

	vector<vector<string>> strings;

	char* file_string =  "C:\\Users\\Hontik\\Desktop\\Bayes\\dataADClosePercDiff"; // "C:\\Users\\Hontik\\Desktop\\Bayes\\ar_normal_single"; // 

	char dfstring[80];
	strcpy(dfstring,file_string);
	strcat(dfstring,".txt");
	
	
	mat data_matrix;
	ifstream myfile(dfstring);
	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 = 50;window_size < 51;window_size++)
		{
			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 fstring[80];					
	strcpy(fstring,file_string);
	strcat(fstring,"lognc.txt");
	//strcat(fstring,"preds.txt");

	for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols() 
	{ 	//pocet stlpcov data_matrix je pocet casovych krokov
		vec cur_res_lognc;
		// vec preds;
		vector<string> nazvy;
		for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();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 (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu?

			if((*model_ref)->my_rarx!=NULL)	//vklada normalizacnz faktor do cur_res_lognc
			{
				cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_rarx->posterior->_ll());				
			}
			else
			{
				cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_arx->_ll());
			}

			pair<vec,vec> predictions = (*model_ref)->predict(20,time,&LapRNG);

			cout << predictions.first << endl << predictions.second << endl;

			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(fstring,ios::app);
			myfilew << avg_prediction << ",";
			myfilew.close();
			*/

			//preds.ins(0,data_matrix.get(0,time+1));
		}

		/*
		myfilew.open(fstring,ios::app);
		myfilew << data_matrix.get(0,time+1) << endl;
		myfilew.close();
		*/

		result_lognc.ins_col(result_lognc.cols(),cur_res_lognc);
		// result_preds.ins_col(result_preds.cols(),preds);

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

		myfilew.open(fstring,ios::app);
		
		// myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
		
		if(time == max_model_order)
		{ 
			for(int i = 0;i<cur_res_lognc.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 << "1";
				}
				else
				{
					myfilew << "0";
				}
					
				if(models[i]->has_constant)
				{
					myfilew << "1";
				}
				else
				{
					myfilew << "0";
				}

				myfilew << ",";
			}

			myfilew << endl;
		}
		
		for(int i = 0;i<cur_res_lognc.size();i++)
		{
			myfilew << cur_res_lognc[i] << ' ';//zmenil som ciarku ze medzeru
		}

		myfilew << endl;				
		
		myfilew.close();
		
}
	


	// EXPERIMENT: One step ahead price prediction. Comparison of classical and robust model using optimal trading
    //             with maximization of logarithm of one-step ahead wealth.

	
		
		/*
		cout << "One experiment finished." << endl;

		ofstream myfile;
		myfile.open("c:\\robust_ar1.txt",ios::app);
		myfile << endl;
		myfile.close();

		myfile.open("c:\\robust_ar2.txt",ios::app);
		myfile << endl;
		myfile.close();*/
	

	//emlig* emlig1 = new emlig(emlig_size);

	//emlig1->step_me(0);
	//emlig* emlig2 = new emlig(emlig_size);
	
	/*
	emlig1->set_correction_factors(4);

	for(int j = 0;j<emlig1->correction_factors.size();j++)
	{
		for(set<my_ivec>::iterator vec_ref = emlig1->correction_factors[j].begin();vec_ref!=emlig1->correction_factors[j].end();vec_ref++)
		{
			cout << j << "    ";
			
			for(int i=0;i<(*vec_ref).size();i++)
			{
				cout << (*vec_ref)[i];
			}

			cout << endl;
		}
	}*/
	
	/*
    vec condition5 = "1.0 1.0 1.01";//"-0.3 1.7 1.5";

	emlig1->add_condition(condition5);
	//emlig1->step_me(0);


	vec condition1a = "-1.0 1.02 0.5";
	//vec condition1b = "1.0 1.0 1.01";
	emlig1->add_condition(condition1a);
	//emlig2->add_condition(condition1b);

	vec condition2a = "-0.3 1.7 1.5";
	//vec condition2b = "-1.0 1.0 1.0";
	emlig1->add_condition(condition2a);
	//emlig2->add_condition(condition2b);

	vec condition3a = "0.5 -1.01 1.0";
	//vec condition3b = "0.5 -1.01 1.0";

	emlig1->add_condition(condition3a);
	//emlig2->add_condition(condition3b);	

	vec condition4a = "-0.5 -1.0 1.0";
	//vec condition4b = "-0.5 -1.0 1.0";	

	emlig1->add_condition(condition4a);
	//cout << "************************************************" << endl;
	//emlig2->add_condition(condition4b);
	//cout << "************************************************" << endl;
	
	//cout << emlig1->minimal_vertex->get_coordinates();
	
	//emlig1->remove_condition(condition3a);
	//emlig1->step_me(0);
	//emlig1->remove_condition(condition2a);
	//emlig1->remove_condition(condition1a);
	//emlig1->remove_condition(condition5);
	

	//emlig1->step_me(0);
	//emlig2->step_me(0);
	

	// DA SE POUZIT PRO VYPIS DO SOUBORU
	// emlig1->step_me(0);

	//emlig1->remove_condition(condition1);
	
	

	
	
	/*
	for(int i = 0;i<100;i++)
	{
		cout << endl << "Step:" << i << endl;		

		double condition[emlig_size+1];		

		for(int k = 0;k<=emlig_size;k++)
		{
			condition[k] = (rand()-RAND_MAX/2)/1000.0;		
		}
			

		vec* condition_vec = new vec(condition,emlig_size+1);
		emlig1->add_condition(*condition_vec);

		/*
		for(polyhedron* toprow_ref = emlig1->statistic.rows[emlig_size]; toprow_ref != emlig1->statistic.end_poly; toprow_ref = toprow_ref->next_poly)
		{
			cout << ((toprow*)toprow_ref)->probability << endl;
		}
		*/
		/*
		cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
	
		/*
		if(i-emlig1->number_of_parameters >= 0)
		{
			pause(30);
		}
		*/

		// emlig1->step_me(i);
		
		/*
		vector<int> sizevector;
		for(int s = 0;s<=emlig1->number_of_parameters;s++)
		{
			sizevector.push_back(emlig1->statistic_rowsize(s));
		}
		*/
	//}
    


	
	/*
	emlig1->step_me(1);

	vec condition = "2.0 0.0 1.0";	

	emlig1->add_condition(condition);

	vector<int> sizevector;
	for(int s = 0;s<=emlig1->number_of_parameters;s++)
	{
		sizevector.push_back(emlig1->statistic_rowsize(s));
	}

	emlig1->step_me(2);

	condition = "2.0 1.0 0.0";

	emlig1->add_condition(condition);

	sizevector.clear();
	for(int s = 0;s<=emlig1->number_of_parameters;s++)
	{
		sizevector.push_back(emlig1->statistic_rowsize(s));
	}
	*/

	return 0;
}