#include <stdio.h>
#include <gsl/gsl_multifit.h>
#include <assert.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>

/** thanks to Jolkdarr, le samedi 18 septembre 2004 à 12:36:37 for 
 * http://www.commentcamarche.net/forum/affich-1009440-lire-une-matrice-dans-un-fichier-texte-en-c**/

using namespace std;


class ligne_matrice : public string
{
};

inline
istream& operator>>(istream& In, ligne_matrice& L)
{
	return getline(In, L);
}




/**    
*/

/**
 *
 *
 *
 gsll_test.m ===================================

#!/usr/bin/octave 

x=randn(50,5);
y=1+x(:,1)+randn;

for i=1:size(x,1),
printf("%4.4f ",y(i));
for j=1:size(x,2),
printf("%d:%4.4f ",j,x(i,j));
end
printf("\n");
end
==================================================
 *
 *** compiling if gsl is installed properly *****************************
   g++ -o gsll gsll.c -lgsl -lgslcblas -lm 
   g++ -o gslp gslp.c -lgsl -lgslcblas -lm 
 *** compiling if gsl is not installed properly (adapt the path to the gsl) ************************
 *** best solution: request gsl on your installation (gsl is just great) ***************************
   g++ -I./gsl-1.9/ -o gsll gsll.c -L./gsl-1.9/.libs -lgsl -lm -Lgsl-1.9/cblas/.libs/ -lgslcblas
   g++ -I./gsl-1.9/ -o gslp gslp.c -L./gsl-1.9/.libs -lgsl -lm -Lgsl-1.9/cblas/.libs/ -lgslcblas
 
    #### you can create the learning set as follows: ./gsll_test.m > learningset 
    ####           only if you have octave --- otherwise create a learning set manually
    or use the learningset at the end of this file (just copy-paste)

    #### learning
   ./gsll learningset model
    #### testing
   ./gslp learningset model outputvalues ### here the learning set is equal to the test set
    #### checking
    paste outputvalues learningset | awk ' { print $1 , " ", $2 }'      #### prints the predicted and the real values
 *
 *
 *
 *
 *  ***/
	int
main (int argc,char*argv[])
{
	char inputFile[5000];
	char outputFile[5000];
	int n;
	int d;
	if (argc<3) {printf("I need two arguments: gsll learningset model \n");exit(-1);}
	sscanf(argv[1],"%s ",&inputFile);
	sscanf(argv[2],"%s ",&outputFile);

/************** reading matrix in stdin ***********************/


	ifstream FichierMatrice(inputFile,ios::in);
	if (!FichierMatrice)
	{
		cerr << "Ouverture du fichier impossible" << endl;
		exit(1);
	}

	double val;double target;char c;int idx;
	unsigned int Lignes = 0;
	unsigned int Colonnes = 0;
	
	std::vector<std::vector<double> > matrix;
	std::vector<double> theTargets;

	// lecture de chaque ligne :
	ligne_matrice Ligne;
	while (FichierMatrice >> Ligne)
	{
		std::vector<double> line;
		Lignes++;
		Colonnes = 0;
		istringstream S(Ligne, ios_base::in);
		S >> target;
		theTargets.push_back(target);
		//cout << "target:" << target << "             ";
		// lecture de chaque element :
 		while (S >> idx)
 		{
 			Colonnes++;
   			//cout <<  "(" << idx << ")  ";
			assert(Colonnes==idx);
			S >> c;
			//cout << "#" << c << "#";
			S >> val;
   			//cout <<  " " << val << "   ";
			line.push_back(val);
 		}
 		matrix.push_back(line);
 	//	cout << endl << "Ligne " << Lignes << " contient " << Colonnes << " colonnes." << endl;		
	}

	n= matrix.size();
	d= matrix[0].size();
	fprintf(stderr,"%d lines, %d columns\n",n,d); 



/************** matrix read ***********************************/

	//GSL::MultiFit::Workspace.alloc(n, p);
	gsl_matrix* x = gsl_matrix_alloc(n,d+1);
	gsl_matrix* cov = gsl_matrix_alloc(d+1,d+1);
	gsl_vector* y=gsl_vector_alloc(n);
	gsl_vector* w=gsl_vector_alloc(d+1);
	fprintf(stderr,"GSL allocations done\n");
	double *chisq=new double(1);
	fprintf(stderr,"double allocation done\n");
	gsl_multifit_linear_workspace *work= gsl_multifit_linear_alloc(n,d+1);
	fprintf(stderr,"main alloc done\n");

	for (int i=0;i<n;i++)
	{
	gsl_vector_set(y,i,theTargets[i]);
	for (int j=0;j<d;j++) gsl_matrix_set(x,i,j,matrix[i][j]);
	gsl_matrix_set(x,i,d,1.);
	}
	fprintf(stderr,"reading done\n");

	gsl_multifit_linear(x,y,w,cov,chisq,work);
	fprintf(stderr,"main job done\n");

	FILE * outputFileH = fopen(outputFile,"w");
	for (int j=0;j<d+1;j++)
	{
		double prout = gsl_vector_get(w,j);
		fprintf(outputFileH,"%g ",prout);
	}
	fprintf(outputFileH,"\n");
	fclose(outputFileH);

	/*
	   int i, n = 4;
	   double x[4] = { 1970, 1980, 1990, 2000 };
	   double y[4] = {   12,   11,   14,   13 };
	   double w[4] = {  0.1,  0.2,  0.3,  0.4 };

	   double c0, c1, cov00, cov01, cov11, chisq;

	   gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
	   &c0, &c1, &cov00, &cov01, &cov11, 
	   &chisq);

	   printf("# best fit: Y = %g + %g X\n", c0, c1);
	   printf("# covariance matrix:\n");
	   printf("# [ %g, %g\n#   %g, %g]\n", 
	   cov00, cov01, cov01, cov11);
	   printf("# chisq = %g\n", chisq);

	   for (i = 0; i < n; i++)
	   printf("data: %g %g %g\n", 
	   x[i], y[i], 1/sqrt(w[i]));

	   printf("\n");

	   for (i = -30; i < 130; i++)
	   {
	   double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
	   double yf, yf_err;

	   gsl_fit_linear_est (xf, 
	   c0, c1, 
	   cov00, cov01, cov11, 
	   &yf, &yf_err);

	   printf("fit: %g %g\n", xf, yf);
	   printf("hi : %g %g\n", xf, yf + yf_err);
	   printf("lo : %g %g\n", xf, yf - yf_err);
	   }
	   return 0;
	   */
}







/****************** example of learning set *************
-14.2881 1:-2.1539 2:-0.8064 3:0.3599 4:1.5678 5:0.5909 
-10.7168 1:-1.6868 2:0.0412 3:0.2989 4:0.5506 5:1.6265 
-1.5478 1:-0.3365 2:1.2474 3:-0.0712 4:-1.5367 5:-0.0756 
8.5313 1:0.9037 2:-1.8271 3:1.1505 4:-0.6216 5:-0.7890 
2.0340 1:0.0534 2:-0.4400 3:1.0541 4:1.1224 5:0.6182 
3.6371 1:0.3744 2:1.4763 3:-0.8212 4:-0.3231 5:-0.4993 
-2.3663 1:-0.3824 2:0.8704 3:-1.2390 4:-0.7142 5:-0.7766 
1.9105 1:0.1738 2:-1.8403 3:0.7640 4:0.1848 5:1.0354 
3.9771 1:0.4907 2:-0.5790 3:-0.6055 4:-0.7748 5:-0.3977 
3.7776 1:0.3897 2:-0.0589 3:-0.4941 4:0.0731 5:0.5806 
-1.5231 1:-0.3566 2:0.3767 3:-0.2675 4:0.3015 5:-0.1446 
2.8262 1:0.2873 2:-0.6907 3:-0.9488 4:0.8091 5:-1.3417 
-8.4526 1:-1.2431 2:1.0786 3:0.2947 4:-0.6611 5:0.8423 
2.9257 1:0.2429 2:-1.6150 3:0.0550 4:-0.2680 5:0.0342 
5.0362 1:0.7817 2:0.2338 3:-0.8249 4:0.1973 5:0.4454 
5.3026 1:0.6739 2:0.7902 3:0.8861 4:-0.2399 5:-0.9348 
-2.7011 1:-0.4794 2:-0.5980 3:-0.7799 4:1.5365 5:-0.7841 
-10.1263 1:-1.5086 2:0.9077 3:0.0606 4:0.7004 5:-0.1054 
2.3707 1:0.1921 2:0.1228 3:1.1855 4:0.2055 5:0.2334 
-7.3536 1:-1.2079 2:0.0329 3:0.6672 4:-0.3958 5:0.5004 
-9.7055 1:-1.5694 2:2.1112 3:0.1536 4:-0.9531 5:0.1273 
-7.9449 1:-1.2742 2:-0.1371 3:-0.3349 4:0.3963 5:-1.5565 
10.7669 1:1.4268 2:-0.4322 3:0.0162 4:1.1060 5:-0.6829 
2.6779 1:0.1457 2:1.1740 3:0.5556 4:0.4947 5:-0.0532 
-6.4095 1:-1.3206 2:-0.8965 3:1.3677 4:-1.4986 5:-0.4771 
2.2105 1:0.1829 2:0.1457 3:-0.8484 4:0.6594 5:-0.1846 
-6.7094 1:-1.2310 2:1.6252 3:1.0419 4:1.1846 5:1.3666 
9.1894 1:1.0182 2:1.5991 3:0.4533 4:0.0958 5:0.2646 
2.0337 1:0.3371 2:1.4003 3:-0.9491 4:-0.7593 5:0.2459 
-7.2570 1:-1.0910 2:-0.2912 3:1.2008 4:0.2384 5:0.7812 
18.1661 1:2.4201 2:-2.2921 3:0.9369 4:-0.0146 5:-0.4881 
3.7458 1:0.5748 2:1.1791 3:-0.2726 4:-0.3721 5:1.9824 
0.8661 1:-0.0993 2:1.0121 3:0.4311 4:0.2627 5:0.9958 
4.7447 1:0.7544 2:-0.0165 3:1.5366 4:0.9906 5:-0.3432 
4.3992 1:0.4946 2:0.6354 3:0.1391 4:1.9843 5:1.4446 
-0.0941 1:-0.1147 2:0.9284 3:0.4452 4:1.1170 5:0.5737 
-4.0578 1:-0.7570 2:0.0238 3:-1.4696 4:1.4182 5:0.9135 
5.3070 1:0.6835 2:0.2403 3:-0.7808 4:-0.4376 5:-1.9839 
6.6623 1:0.9085 2:0.5257 3:1.0482 4:0.2890 5:-0.1438 
4.4416 1:0.7135 2:0.0165 3:-1.1667 4:-0.7540 5:0.0439 
-3.6725 1:-0.7476 2:0.9952 3:-0.9166 4:-1.0452 5:-1.3816 
-0.9815 1:-0.4563 2:-0.4585 3:1.8102 4:-0.5970 5:0.6710 
-2.9490 1:-0.5859 2:-0.6639 3:-0.7555 4:0.8469 5:-0.3758 
3.9295 1:0.6190 2:-0.5269 3:-0.5739 4:-0.9378 5:-1.3872 
12.2029 1:1.5079 2:0.1130 3:-1.3663 4:-0.3483 5:0.3997 
7.1899 1:1.0230 2:-0.7200 3:-1.7784 4:1.1260 5:-0.1766 
9.9616 1:1.0496 2:0.8701 3:0.4646 4:-0.1699 5:-1.4635 
-17.6241 1:-2.4989 2:1.5232 3:1.5566 4:0.4926 5:0.2047 
-7.1149 1:-1.0837 2:-0.6059 3:0.1088 4:-0.2093 5:-0.4709 
12.3284 1:1.6331 2:-0.8860 3:0.7226 4:1.6387 5:-1.1022 
***************************/