#include "optim.h"

PVecteurRn allouerVecteurRn(int dimension) {
	if(dimension <= 0)
		return (PVecteurRn) NULL;
	else {
		PVecteurRn pv;
		int i;

		pv = (PVecteurRn) malloc(sizeof(SVecteurRn));
		pv->dimension = dimension;
		pv->pvecteur = (double *) malloc(dimension * sizeof(double));
		for (i = 0; i < dimension;i++)
			pv->pvecteur[i] = 0.0;
		return pv;
	}
}

PMatriceRn allouerMatriceRn(int dimension) {
	if(dimension <= 0)
		return (PMatriceRn) NULL;
	else {
		PMatriceRn pm;
		int i,j;

		pm = (PMatriceRn) malloc(sizeof(SMatriceRn));
		pm->dimension = dimension;
		pm->pmatrice = (double **) malloc(dimension * sizeof(double *));
		for (i = 0; i < dimension;i++) {
			pm->pmatrice[i] = (double *) malloc(dimension * sizeof(double));
			for(j = 0; j < dimension; j++)
				pm->pmatrice[i][j] = 0.0;
		}
		return pm;
	}
}

PMatriceRn allouerMatriceIdentiteRn(int dimension) {
	if(dimension <= 0)
		return (PMatriceRn) NULL;
	else {
		PMatriceRn pm;
		int i,j;

		pm = (PMatriceRn) malloc(sizeof(SMatriceRn));
		pm->dimension = dimension;
		pm->pmatrice = (double **) malloc(dimension * sizeof(double *));
		for (i = 0; i < dimension;i++) {
			pm->pmatrice[i] = (double *) malloc(dimension * sizeof(double));
			for(j = 0; j < dimension; j++)
                if(i==j){
                    pm->pmatrice[i][j] = 1.0;
                }else{
                    pm->pmatrice[i][j] = 0.0;
                }
		}
		return pm;
	}
}

PVecteurRn copierVecteurRn(PVecteurRn pvect) {
	int dimension = pvect->dimension;
	PVecteurRn pres = allouerVecteurRn(dimension);
	int i;

	for (i = 0; i < dimension;i++)
		pres->pvecteur[i] = pvect->pvecteur[i];
	return pres;
}

PMatriceRn copierMatriceRn(PMatriceRn pmat) {
	int dimension = pmat->dimension;
	PMatriceRn pres = allouerMatriceRn(dimension);
	int i,j;

	for (i = 0; i < dimension;i++)
		for(j = 0; j < dimension; j++)
			pres->pmatrice[i][j] = pmat->pmatrice[i][j];
	return pres;
}


void libererVecteurRn(PVecteurRn pvect) {
	free(pvect->pvecteur);
	free(pvect);
}

void libererMatriceRn(PMatriceRn pmat) {
	int dimension = pmat->dimension;
	int i;

	for(i = 0; i < dimension;i++)
		free(pmat->pmatrice[i]);
	free(pmat->pmatrice);
	free(pmat);
}

PVecteurRn addVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	PVecteurRn pvect3;

	if(pvect1->dimension != pvect2->dimension) {
		return (PVecteurRn)NULL;
	}
	else {
		int i;
		int dimension = pvect1->dimension;

		pvect3 = allouerVecteurRn(dimension);
		for(i = 0; i < dimension;i++) {
			pvect3->pvecteur[i] = pvect1->pvecteur[i] + pvect2->pvecteur[i];
		}
		return pvect3;
	}
}

PMatriceRn addMatriceRn(PMatriceRn pmat1, PMatriceRn pmat2) {
	PMatriceRn pmat3;

	if(pmat1->dimension != pmat2->dimension) {
		return (PMatriceRn)NULL;
	}
	else {
		int i,j;
		int dimension = pmat1->dimension;

		pmat3 = allouerMatriceRn(dimension);
		pmat3->dimension = dimension;
		for(i = 0; i < dimension;i++) {
			for(j = 0; j < dimension;j++) {
				pmat3->pmatrice[i][j] = pmat1->pmatrice[i][j] + pmat2->pmatrice[i][j];
			}
		}
		return pmat3;
	}
}

PMatriceRn mulMatricesRn(PMatriceRn pmat1, PMatriceRn pmat2) {
	PMatriceRn pmat3;

	if(pmat1->dimension != pmat2->dimension) {
		return (PMatriceRn)NULL;
	}
	else {
		int i,j,k;
		int dimension = pmat1->dimension;

		pmat3 = allouerMatriceRn(dimension);
		pmat3->dimension = dimension;
		for(i = 0; i < dimension;i++) {
			for(j = 0; j < dimension;j++) {
				double res;

				res = 0.0;
				for(k = 0; k < dimension;k++)
					res += pmat1->pmatrice[i][k] + pmat2->pmatrice[k][j];
				pmat3->pmatrice[i][j] = res;
			}
		}
		return pmat3;
	}
}


PVecteurRn difVecteurRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	PVecteurRn pvect3;

	if(pvect1->dimension != pvect2->dimension) {
		return (PVecteurRn)NULL;
	}
	else {
		int i;
		int dimension = pvect1->dimension;

		pvect3 = allouerVecteurRn(dimension);
		pvect3->dimension = dimension;
		pvect3->pvecteur =(double *) malloc(dimension * sizeof(double));
		for(i = 0; i < dimension;i++) {
			pvect3->pvecteur[i] = pvect1->pvecteur[i] - pvect2->pvecteur[i];
		}
		return pvect3;
	}
}

double produitScalaireRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	if(pvect1->dimension != pvect1->dimension)
		return 0.0;
	else {
		double res = 0.0;
		int dimension = pvect1->dimension;
		int i;

		for(i = 0; i < dimension; i++)
			res += pvect1->pvecteur[i] * pvect2->pvecteur[i];
		return res;
	}
}

PMatriceRn produitVecteursRn(PVecteurRn pvect1, PVecteurRn pvect2) {
	if(pvect1->dimension != pvect1->dimension)
		return (PMatriceRn)NULL;
	else {
		int dimension = pvect1->dimension;
		PMatriceRn pm = allouerMatriceRn(dimension);

		int i,j;

		for(i = 0; i < dimension; i++)
			for(j = 0; j < dimension; j++)
				pm->pmatrice[i][j] = pvect1->pvecteur[i] * pvect2->pvecteur[j];
		return pm;
	}
}

double normeCarree(PVecteurRn pvect) {
	int i;
	int dimension = pvect->dimension;
	double res = 0.0;
	double a;

	for(i = 0; i < dimension; i++) {
		a = pvect->pvecteur[i];
		res += a * a;
	}
	return res;
}

double norme(PVecteurRn pvect) {
	return sqrt(normeCarree(pvect));
}


PVecteurRn mulScalVecteurRn(double coeff,PVecteurRn pvect) {
	PVecteurRn pvect2;
	int i;
	int dimension = pvect->dimension;

	pvect2 = allouerVecteurRn(dimension);
	for(i = 0; i < dimension;i++) {
		pvect2->pvecteur[i] = coeff *  pvect->pvecteur[i];
	}
	return pvect2;
}

PMatriceRn mulScalMatriceRn(double coeff,PMatriceRn pmat) {
	PMatriceRn pmat2;
	int i,j;
	int dimension = pmat->dimension;

	pmat2 = allouerMatriceRn(dimension);
	for(i = 0; i < dimension;i++) {
		for(j = 0; j < dimension; j++) {
			pmat2->pmatrice[i][j] = coeff *  pmat->pmatrice[i][j];
		}
	}
	return pmat2;
}


void affVecteurRn(PVecteurRn pvect) {
	int i;
	int dimension = pvect->dimension;

	printf("(");
	for(i = 0; i < dimension;i++) {
		printf("%f ",pvect->pvecteur[i]);
		if(i != (dimension-1))
			printf(",");
	}
	printf(")\n");
}

void affMatriceRn(PMatriceRn pmat) {
	int i,j;
	int dimension = pmat->dimension;

	for(i = 0; i < dimension;i++) {
		for(j=0; j < dimension; j++)
			printf("%f ",pmat->pmatrice[i][j]);
		printf("\n");
	}
}

PVecteurRn produitMatriceVecteurRn(PMatriceRn pmat, PVecteurRn pvect) {
	if(pmat->dimension != pvect->dimension)
		return (PVecteurRn) NULL;
	else {
		int dimension = pmat->dimension;
		PVecteurRn pvectRes = allouerVecteurRn(dimension);
		double res;
		int i,j;

		for(i = 0; i < dimension; i++) {
			res = 0.0;

			for(j = 0; j < dimension ; j++)
				res += pmat->pmatrice[i][j] * pvect->pvecteur[j];
			pvectRes->pvecteur[i] = res;
		}
		return pvectRes;
	}
}


int main() {
	PVecteurRn px0 = allouerVecteurRn(2);
	PVecteurRn pvect;

	px0->pvecteur[0] = 7.0;
	px0->pvecteur[1] = 1.5;
	pvect = gradientPasOptimal(px0,10,.000001,100000);
	affVecteurRn(pvect);
	libererVecteurRn(pvect);
	pvect = gradientPasFixe(px0,0.001,.000001,100000);
	affVecteurRn(pvect);
	libererVecteurRn(pvect);
	return 0;
}