/** @(#) numericalap/forwardgradient.cpp */

// #include "forwardgradient.h"
/** @(#) forwardgradient.h */
#ifndef FORWARDGRADIENT_H
#define FORWARDGRADIENT_H
//----------------------------------------------------------------------------
#include "../functionnd/functionnd_1d.h"
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/vectordouble.h"
#include "../mathpos/mathpos.h"
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include <stdexcept>	// runtime_exception
using namespace std;

void forwardGradient(int n, VectorDouble &xc, double fc,
		FunctionND_1D &theFunction, VectorDouble &Sx,
		VectorDouble &g, double eta);

void checkForwardGradient(int n, VectorDouble &xc, double fc,
		FunctionND_1D &theFunction, VectorDouble &Sx,
		VectorDouble &g, double eta)
                throw(invalid_argument);
/**
* Calcula una aproximacion hacia adelante hacia el gradiente
* grad(f(xc)) usando solo los valores de f(xc).
* El algoritmo es identico al A5.4.1 (FDJAC) excepto que
* los parametros fc (escalar)y g (vector) del algoritmo A5.6.3 son remplazados
* por los parametros Fc (vector) y J (Jacobiano) en el algoritmo A5.4.1
* y la linea 2.6 del algoritmo A.6.3
* son remplazados por el correspondiente ciclo del algoritmo A5.4.1.
* @param n Entero. Natural. Dimension de las matrices y vectores.
* @param xc Vector con punto de prueba inicial.
* @param fc Evaluacion en xc (equivale a f(xc))
* @param theFunction Nombre de la funcion. F:R^n->R. FN
* @param Sx Escala de xc.
* @param g return Aproximacion al grad(f(xc)).
* @param eta Precision de la maquina. Omision: mathpos::MACHEPS.
* @author Omar Posada Villarreal, Daniel Alba Cuellar
* @version 1.0, 04/03/2002
* @see Dennis. Pagina 320
* Algoritmo A5.6.3 (FDGRAD).
* Ejemplo 4.1.3 pag 72
*/
void forwardGradient(int n, VectorDouble &xc, double fc,
		FunctionND_1D &theFunction, VectorDouble &Sx,
		VectorDouble &g, double eta = MACHEPS) {
        //cout << "\n{ forwardGradient";
	// { logica
        xc.setLogic(1);
        Sx.setLogic(1);
	g.setLogic(1);
	// } logica

        int i, j;
        double stepsizej;	// taman~o de paso
        double tempj;
        double absXc;
        double invSx;
        double fj;

        // ALGORITMO
	double sqrteta = sqrt(eta);		// 1.  tolerancia
        for (j = 1; j <= n; ++j) {		// 2.

        	// 2.1 Calcula columna j de J
		absXc = fabs(xc[j]);
		invSx = 1.0 / Sx[j];	// div / 0

                // Se puede cambiar esta regla para diferente taman~o de paso.
                stepsizej =  sign( xc[j] ) * sqrteta * max(absXc, invSx);
                tempj = xc[j];  		// 2.2

                // Reduce ligeramente los errores de prec finita
                xc[j] += stepsizej;     	// 2.3
		stepsizej = xc[j] - tempj;	// 2.4

                // fj = f(xc + stepsizej * ej)
                fj = theFunction.evaluate(xc);	// 2.5

                // { Diferente de A5.4.1 FDJAC
                g[j] = (fj - fc) / stepsizej;	// 2.6
                // } Diferente de A5.4.1 FDJAC

                // Restaura xc[j]
                xc[j] = tempj;	// 2.7
        }
        //cout << "\n  } forwardGradient";
}

/** Checa parametros: dimensiones, rangos, divisiones entre cero */
void checkForwardGradient(int n, VectorDouble &xc, double fc,
		FunctionND_1D &theFunction, VectorDouble &Sx,
		VectorDouble &g, double eta = MACHEPS)
                throw(invalid_argument) {
	// Condiciones invalidas
        if ( (n < 1)
        	|| (xc.getSize() != n)
        	|| (Sx.getSize() != n)
        	|| (g.getSize() != n)
                || (eta < 0.0)
        	) {
        	throw invalid_argument("\nforwardGradient: arg invalidos.");
        }
        // Prevee division entre cero
        if (Sx.hasAlmostZeros()) {
        	throw invalid_argument("forwardGradient: Division entre cero.");
        }

	forwardGradient(n, xc, fc, theFunction, Sx, g, eta);
}

//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------