/** @(#) numericalap/finitenewtondogleg.cpp */
//#include "finitenewtondogleg.h"
/** @(#) finitenewtondogleg.h */
#ifndef FINITENEWTONDOGLEG_H
#define FINITENEWTONDOGLEG_H
//----------------------------------------------------------------------------
#include "../adjustedrecipes/qrdcmp.cpp"
#include "../adjustedrecipes/qrsolv.cpp"
#include "../functionnd/functionnd.h"
#include "../functionnd/functionnd_1d.h"
#include "../matrixdouble/squaredmatrix.h"
#include "../mathpos/mathpos.h" // sign, MACHEPS
#include "../utilpos/interfacepos.h"
#include "factorizecholesky.cpp"
#include "doglegdriver.cpp"
//#include "finitejacobi.cpp";
#include <iostream> // salida
#include <fstream>
#include <algorithm> // max
#include <cmath>
#include <iostream>
#include <stdexcept> // runtime_exception
using namespace std;
void checkFiniteNewtonDogLeg(int n, VectorDouble &xc, double &fc,
FunctionND_1D &theFunction, VectorDouble &Sx,
int iterations,
double maxstep, double steptol,
double &delta,
double eta = MACHEPS,
int logType = SHORT_LOG, int waitType = NO_LOG,
string pathFile = "finiteNewtonDogleg.txt")
throw (invalid_argument, runtime_error);
void finiteNewtonDogLeg(int n, VectorDouble &xc, double &fc,
FunctionND_1D &theFunction, VectorDouble &Sx,
int iterations,
double maxstep, double steptol,
double &delta,
double eta = MACHEPS,
int logType = SHORT_LOG, int waitType = NO_WAIT,
string pathFile = "finiteNewtonDogleg.txt")
throw (runtime_error);
class FunctionG : public FunctionND {
public:
// No se checa las dimensiones de xc y Fc.
void evaluate(VectorDouble &xc, VectorDouble &Fc) {
xc.setLogic(1);
Fc.setLogic(1);
// Funcion del usuario
Fc[1] = 4.0*pow(xc[1], 3.0) + 2.0*xc[1];
Fc[2] = 2.0*xc[2];
}
}; // No olvidar ";"
void hessianEvaluate(VectorDouble &xc, SquaredMatrix &H) {
xc.setLogic(1);
H.setLogic(1);
// Funcion del usuario
H[1][1] = 12.0 * pow(xc[1], 2.0) + 2.0;
H[1][2] = 0.0;
H[2][1] = 0.0;
H[2][2] = 2.0;
}
/** Aproximacion por diferencias finitas por el metodo de quasi Newton
* Para la decision del paso elige entre elpaso de Newton o el DogLeg.
* No se checan dimensiones. La dimension n se toma de la dimension de xc.
* Logica uno.
* @author Daniel Alba Cuellar, Omar Posada Villarreal
* @version 1.0, 22/04/2002
* @param n Entero. Natural. Dimension de las matrices y vectores.
* @param xc IN/OUT
* IN: Vector con punto de prueba inicial.
* OUT: Solucion
* @param fc IN/OUT
* IN: Vector con evaluacion punto de prueba inicial.
* OUT: Evaluacion de la solucion
* @param theFunction Nombre de la funcion. F:R^n->R^n. FVEC.
* @param Sx Escala. Valores usuales.
* @param maxstep Maximo taman~o de paso. Dado por QuasiNewton
* @param steptol Precision del paso
* @param delta in/out Radio de la region de confianza.
* @param xp return x+
* @param fp return f+
* @param eta Precision de la maquina. Omision: mathpos::MACHEPS.
* @param deb debug. Mostrar resultados.
* 0: No mostrar.
* 1: Mostrar resumen, version corta.
* 2: Mostrar pasos intermedios, version larga.
* @param logType
* @see pages 314-315
*/
/*
* @param H return. Matriz cuadrada Hessiana.
* @param fc return Vector sobreescrito Entrada/Salida.
*/
void finiteNewtonDogLeg(int n, VectorDouble &xc, double &fc,
FunctionND_1D &theFunction, VectorDouble &Sx,
int iterations,
double maxstep, double steptol,
double &delta,
double eta,
int logType, int waitType,
string pathFile)
throw (runtime_error) {
// QUITAR MAX_LOG Y OTROS, CAMBIAR LA FIRMA DE LAS FUNCIONES
// En depuracion usar checkMethod, al final del archivo
LLOG << "\n{ finiteNewton";
// log file
ofstream outFile(pathFile.data(), ios::out);
if ( !outFile ) { // !sobrecargado
throw runtime_error(
"finiteNewtonDogLeg: Error de apertura de archivo.");
}
// { posada
xc.setLogic(1);
Sx.setLogic(1);
//H.setLogic(1);
// } posada
FunctionG funcGrad;
SquaredMatrix H(n, 1); // Hessiana
SquaredMatrix R(n, 1); // Factoriz QR
SquaredMatrix L(n, 1); // Factoriz Cholesky
VectorDouble c(n, 1); // used by qrdcmp
VectorDouble d(n, 1); // used by qrdcmp
VectorDouble g(n, 1); // grad(f(xk))
VectorDouble xp(n, 1); // x+
VectorDouble sN(n, 1); // en el libro lo llama sk^N
// param de IN/OUT, en Ax=b en IN es b, y en OUT x
//sN Solucion aproximada a H(f(xk)) sk^N = - grad(f(xk))
bool maxtaken; // x+ esta en el radio de la region de confinza
int retcode;
int i, k;
int sing;
double cond;
double mx, num, den, aux; // condicion de paro
double fp; // f+
bool Newttaken;
// Encabezado
SLOG << "Iter\tRadio1\tRadio2\tPaso\tfc\tdim(xc)\txc";
FLOG << "Iter\tRadio1\tRadio2\tPaso\tfc\tdim(xc)\txc";
// iteracion 0
SLOG << "\n0\t" << delta << "\t" << delta << "\tInicial"
<< "\t" << fc << "\t" << n << "\t" << xc;
FLOG << "\n0\t" << delta << "\t" << delta << "\tInicial"
<< "\t" << fc << "\t" << n << "\t" << xc;
fc = theFunction.evaluate(xc); // antes param
// Valor tipico de f typ f
double typF = theFunction.evaluate(Sx);
// newton iteration
for (k=1; k<=iterations; k++) {
funcGrad.evaluate(xc, g);
hessianEvaluate(xc, H);
//// finiteJacobi(n, xc, Fc, theFunction, Sx, J, eta); //aprox por dif fin
//// Fc.multiply(-1.0, negFc);
sing = qrdcmp(H, n, c, d); // aproximacion
// Si es singular, no se puede buscar la condicion.
if (sing) {
throw runtime_error(
"\nfiniteNewtonDogLeg: La matriz es singular.");
} else {
//sN Solucion aproximada a H(f(xk)) sk^N = -grad(f(xk))
R.createRMatrix(H, d);
/* qrdcmp(fQR, n, c, fd, &sing);
createRMatrix(QR, d, R, n);
double cond = condest(R, n);
*/
cond = R.condition();
if (cond > (1.0 / eta)) {
throw runtime_error(
"\nfiniteNewtonDogLeg: Hessiano mal condicionado.");
}
sN = g * -1.0;
qrsolv(H, n, c, d, sN); //sN in/out
// Obtener triang inf de hessiana
if ( !factorizeCholesky(n, H, L)) {
throw runtime_error(
"\nfiniteNewtonDogLeg: No se pudo factorizar H.");
}
// delta es in/out
SLOG << "\n" << k << "\t" << delta; // IN
FLOG << "\n" << k << "\t" << delta; // IN
// en dogdriver se actualiza x
doubleDogLegDriver(n, xc, fc,
theFunction, g, L,
sN, Sx,
maxstep, steptol,
delta,
retcode, xp, fp, maxtaken, Newttaken);
SLOG << "\t" << delta; // OUT
FLOG << "\t" << delta; // OUT
if (maxtaken) {
LLOG << "\nx+ esta cerca del"
<< " radio de la region de confinza.";
}
// duda con las evaluaciones de fp y fc
if (Newttaken) {
// Usa el paso de Newton
xc = xc + sN;
fc = theFunction.evaluate(xc);
SLOG << "\tNewton";
FLOG << "\tNewton";
} else {
// Usa el doglegstep
xc = xp;
fc = fp;
SLOG << "\tDogleg";
FLOG << "\tDogleg";
}
/* 0: x+ satisfactorio encontrado
* 1: la rutina fallo para encontrar una x+ satisfactoria,
* suficientemente distinta de xc.
*/
//if (retcode == 0) { // No hacer nada
if (retcode == 1) {
LLOG << "\n } finiteNewtonDogLeg";
return;
}
SLOG << "\t" << fc << "\t" << n << "\t" << xc;
FLOG << "\t" << fc << "\t" << n << "\t" << xc;
// Criterio de paro
// Dennis. pag 160. Ecuacion (7.2.5)
// max (i=[1, n]) | num/den | <= steptol
// num = grad(f(x)i max{|(x+)i|, typ xi}
// den = max{|f(x+)|, typ f}
// typf = f(xc)
den = max( fabs(fc), fabs(typF));
if (isAlmostZero(den)) {
throw runtime_error(
"\nfiniteNewtonDogLeg: Division entre cero.");
}
mx = -1; // mx debe ser positivo
for (i = 1; i <= n; ++i) {
num = g[i] * max( fabs(xp[i]), 1.0/Sx[i] );
aux = fabs(num / den);
if (aux > mx) {
mx = aux;
}
}
LLOG << "\nCondicion de paro. mx < steptol ("
<< mx << " < " << steptol << "): "
<< toStr(mx < steptol);
if (mx < steptol) {
LLOG << "\n } finiteNewtonDogLeg";
return;
}
SWAIT();
} // if 1
}
LLOG << "\n } finiteNewtonDogLeg";
}
/** Checa parametros: dimensiones, rangos. */
void checkFiniteNewtonDogLeg(int n, VectorDouble &xc, double &fc,
FunctionND_1D &theFunction, VectorDouble &Sx,
int iterations,
double maxstep, double steptol,
double &delta,
double eta,
int logType, int waitType,
string pathFile)
throw (invalid_argument, runtime_error) {
// Condiciones invalidas
if ( (n < 1)
|| (xc.getSize() != n)
|| (Sx.getSize() != n)
|| (iterations < 1)
|| (delta < 0.0)
|| (eta < 0.0)
) {
throw invalid_argument("\nfiniteNewtonDogLeg: arg invalidos.");
}
// Prevee division entre cero
if (Sx.hasAlmostZeros()) {
throw invalid_argument(
"finiteNewtonDogLeg: Division entre cero.");
}
finiteNewtonDogLeg(n, xc, fc, theFunction, Sx, iterations,
maxstep, steptol, delta, eta,
logType, waitType, pathFile);
}
//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------