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

// #include "finitehessiangradients.h"
/** @(#) finitehessiangradients.h */
#ifndef FINITEHESSIANFUNCTION_H
#define FINITEHESSIANFUNCTION_H
//----------------------------------------------------------------------------
#include "finitejacobi.cpp"
#include "../functionnd/functionnd.h"
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/vectordouble.h"
#include "../mathpos/mathpos.h"
#include <iomanip>
#include <iostream>
#include <string>
#include <stdexcept>	// runtime_exception
using namespace std;

/**
* Calcula una aproximacion hacia adelante mediante diferencias finitas hacia
* la matriz Hessiana (grad^2(f(xc)) usando valores analiticos
* del gradiente (grad(f(xc)).
* Uso:
*        func.evaluate(xc, g);
*	finiteHessianGradients(n, xc, g, func, Sx, H);
* Solo la diagonal y el triangulo superior de H tienen valores correctos
* despues de regresar del algoritmo, debido a que estas son las unicas
* porciones de H que son referenciadas por el resto del
* sistema de algoritmos. De cualquier forma, le matriz completa H es usada
* como un paso intermedio.
* No necesita cambios, si se desea economizar el almacenaminto.
* @param xc Vector con punto de prueba inicial. De aqui se toma la dimension
* @param grad Funcion gradiente. grad(f) : R^n -> R^n
* @param g grad(f(xc)) Evaluacion del gradiente.
* @param Sx Escala de xc.
* @param H return. Matriz hessiana simetrica ( equiv grad^2(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.1 (FDHESSG).
* Ejemplo 4.1.3  pag 72.
*/
void finiteHessianGradients(int n, VectorDouble &xc, VectorDouble &g,
		FunctionND &grad, VectorDouble &Sx,
		SquaredMatrix &H, double eta = MACHEPS)
                throw (invalid_argument) {
	finiteJacobi(n, xc, g, grad, Sx, H, eta);
	int i, j;
	// H = (H + Ht) / 2
	for (i = 1; i <= (n - 1); ++i) {
		for (j = (i + 1); j <= n; ++j) {
			H[i][j] = (H[i][j] + H[j][i]) / 2.0;
		}
	}
}

/** Checa parametros: dimensiones, rangos, divisiones entre cero */
void checkFiniteHessianGradients(int n, VectorDouble &xc, VectorDouble &g,
		FunctionND &grad, VectorDouble &Sx,
		SquaredMatrix &H, double eta = MACHEPS)
                throw (invalid_argument) {
	// Condiciones invalidas
        if ( (n < 1)
        	|| (xc.getSize() != n)
        	|| (g.getSize() != n)
        	|| (Sx.getSize() != n)
        	|| (H.getSize() != n)
                || (eta < 0.0)
        	) {
        	throw invalid_argument("\nfiniteHessFunc: arg invalidos.");
        }
        // Prevee division entre cero
        if (Sx.hasAlmostZeros()) {
        	throw invalid_argument("finiteHessFunc: Division entre cero.");
        }
	finiteHessianGradients(n, xc, g, grad, Sx, H, eta);
}

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