#include <stdio.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.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);
}




/**    g++ -I./gsl-1.9/ -o gslp gslp.c -L./gsl-1.9/.libs -lgsl -lm -Lgsl-1.9/cblas/.libs/ -lgslcblas
*/

/**
 *
 *
 *
 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
==================================================
 *
 *
 *  
 *
 *  ***/
	int
main (int argc,char*argv[])
{
	char inputFile[5000];
	char outputFile[5000];
	char model[5000];
	int n;
	int d;
	if (argc<4) {printf("I need three arguments: gslp testset model outputvalues\n");exit(-1);}
	sscanf(argv[1],"%s ",&inputFile);
	sscanf(argv[2],"%s ",&model);
	sscanf(argv[3],"%s ",&outputFile);

/************** reading matrix in inputFile ***********************/

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

	//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* yp=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,"matrix reading done\n");


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


/************** reading vector in model ***********************/

	FILE *modelFile = fopen(model,"r");
	for (int i=0;i<d+1;i++)
	{
		double prout;
		fscanf(modelFile,"%lf",&prout);
		gsl_vector_set(w,i,prout);
	}
	fclose(modelFile);

/************** vector read ***********************************/

	gsl_blas_dgemv( CblasNoTrans, 1.0, x, w, 0.0, yp );

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

}