#include "optim.h"

/**
Les coefficients utiles dans la règle d'Armijo
*/
#define TAU 0.01
#define OMEGA1 0.0001

double pasArmijo(PVecteurRn px, PVecteurRn pdirectionDescente, double pas) {
	PVecteurRn pdescente;
	PVecteurRn pxNew;
	double JxOld = J(px);
	double grCarre = normeCarree(pdirectionDescente);

	pdescente = mulScalVecteurRn(pas,pdirectionDescente);
	pxNew = addVecteurRn(px,pdescente);
	//affVecteurRn(px);
	while(J(pxNew) > (JxOld - OMEGA1 * pas * grCarre)) {
		//printf("norme %f  J(pxNew)=%f   JxOld=%f\n",grCarre,J(pxNew), JxOld);
		pas /= 2.0;
		libererVecteurRn(pdescente);
		libererVecteurRn(pxNew);
		pdescente = mulScalVecteurRn(pas,pdirectionDescente);
		pxNew = addVecteurRn(px,pdescente);
		//affVecteurRn(pxNew);
	}
	//printf("\n\n");
	if (pas < 0.001) return 0.001;
	return pas;
}

PVecteurRn gradientPasOptimal(PVecteurRn px0, double pas,double epsilon, int nbMaxIterations) {
    PVecteurRn pxOld;
    PVecteurRn pxNew;
    PVecteurRn pgradient;
	PVecteurRn pdescente;
	PVecteurRn pdirectionDescente;
	PVecteurRn pdiff;
	int nbIterations;
	double ecart;
	double pas1;

	pxOld = copierVecteurRn(px0);

	nbIterations = 0;

	do {
		//printf("iter %d\n",nbIterations +1);
		//printf("x%d ",nbIterations);affVecteurRn(pxOld);
		pgradient = gradientJ(pxOld);
		pdirectionDescente = mulScalVecteurRn(-1.0,pgradient);


		pas1 = pasArmijo(pxOld,pdirectionDescente,pas);
		pdescente = mulScalVecteurRn(pas1,pdirectionDescente);
		if(nbIterations <= 3)
			affVecteurRn(pdescente);
		//printf("descente%d ",nbIterations);affVecteurRn(pdescente);
		//printf("\n\n");
		pxNew = addVecteurRn(pxOld,pdescente);
		pdiff = difVecteurRn(pxNew,pxOld);
		ecart = norme(pdiff);

		libererVecteurRn(pgradient);
		libererVecteurRn(pdirectionDescente);
		libererVecteurRn(pdescente);
		libererVecteurRn(pdiff);
		libererVecteurRn(pxOld);
		pxOld = pxNew;
		//affVecteurRn(pxOld);
		nbIterations++;
	} while (ecart >= epsilon && nbIterations <= nbMaxIterations);
	printf("Pas optimal : nbIter %d ",nbIterations);

	return pxNew;
}