/** @(#) doglegdriver.cpp */
// #include "ap.h"
/** @(#) .h */
/**
* @author Daniel Alba Cuellar
* @author Omar Posada Villarreal
* @version 1.0, /03/2002
*/
#ifndef DOGLEGDRIVER_H
#define DOGLEGDRIVER_H
//----------------------------------------------------------------------------
#include "doglegstep.cpp"
#include "updatetrustregion.cpp"
#include "../functionnd/functionnd_1d.h"
#include "../mathpos/mathpos.h"
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/vectordouble.h"
#include <iomanip>
#include <iostream>
#include <string>
#include <stdexcept> // runtime_exception
using namespace std;
/** Busca una x+ en la curva double dogleg tal que
* f(x+) <= f(xc) + ag^T (x+ - xc)
* (se usa alfa = 10^-4) y un taman~o de paso escalado steplength = delta,
* iniciando con la entrada delta pero incrementando o decrementando
* delta si es necesario. Tambien, produce una region de confianza
* inicial delta para la siguiente iteracion.
* Se llama al algoritmo A6.4.4. y A6.4.5
* A6.4.4. para el taman~o de paso.
* A6.4.5. Aceptar
* Usa logica uno.
* @author Omar Posada Villarreal, Daniel Alba Cuellar
* @param n Taman~o de vectores y matrices. No se checan dimensiones.
* @param xc Vector con punto de prueba inicial. De aqui se toma la dimension
* @param fc Evaluacion de la funcion en xc.
* @param theFunction Nombre de la funcion. F:R^n->R. FN.
* @param g Gradiente
* @param L Triangular inferior
* @param sN = - (L L^T)^-1 g). Neton step.
* @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 retcode return
* 0: x+ satisfactorio encontrado
* 1: la rutina fallo para encontrar una x+ satisfactoria,
* suficientemente distinta de xc.
* @param xp return x+
* @param fp return f+
* @param maxtaken return
* @param Newttaken return Este parametro es solo informativo,
* no vienen en Dennis, se usa para escribir en archivo.
* true: Se usa el paso de Newton
* false: se usa el dogleg step
* @version 1.0, 20/04/2002, sig vers LowerTriangularMatrix &L
* @see Dennis. A6.4.3. DOGDRIVER pagina 335
*/
void doubleDogLegDriver(int n, VectorDouble &xc, double fc,
FunctionND_1D &theFunction, VectorDouble &g, SquaredMatrix &L,
VectorDouble &sN, VectorDouble &Sx,
double maxstep, double steptol,
double &delta,
int &retcode, VectorDouble &xp, double &fp, bool &maxtaken,
bool &Newttaken,
int logType = SHORT_LOG, int waitType = NO_WAIT)
throw (invalid_argument) {
LLOG << "\n{ doubleDogLegDriver";
// didac
LLOG << "\n-----IN-----";
LLOG << "\nDimension. n = " << n;
LLOG << "\nPunto inicial. xc: " << xc;
LLOG << "\nf(xc). Funcion en el punto inicial. fc: " << fc;
LLOG << "\ngrad(f(xc)). g = " << g;
LLOG << "\nTriangular inferior del Hessiano. L:\n" << L;
LLOG << "\nsc^N = -Hc^-1 grad(f(xc)). sN:" << sN;
LLOG << "\nEscala metrica. Sx: " << Sx;
LLOG << "\nMaxima delta. maxstep = " << maxstep;
LLOG << "\nsteptol = " << steptol;
LLOG << "\n-----IN/OUT ANTES-----";
LLOG << "\ndelta c. Radio de la region de confianza. delta = "<< delta;
LWAIT();
// Misma logica
xc.setLogic(1);
g.setLogic(1);
L.setLogic(1);
sN.setLogic(1);
Sx.setLogic(1);
xp.setLogic(1);
// Var de A6.4.3
// Var de A6.4.4
// in
double Newtlen; // init por 3.
// in/out
bool firstdog;
double Cauchylen;
double eta;
VectorDouble sSD(n, 1);
VectorDouble v(n, 1);
// out
VectorDouble s(n, 1);
//bool Newttaken; // Si no se usa param, descomentar
// Var de A6.4.5
// Usado por A6.4.4. VectorDouble &s,
// Usado por A6.4.4. bool Newtakken,
int steptype;
SquaredMatrix H(n, 1); // Dummy, usado para hookstep
VectorDouble xpprev(n, 1);
double fpprev; // ??? init
VectorDouble Dx_sN(n, 1);
// Es una matriz cuadrada, pero se almacenara como un vector a Sx
VectorDouble Dx(Sx);
// Significa una llamada inicial a A6.4.4
steptype = 2; // 0. Usar a dogstep
// A6.4.5.
retcode = 4; // 1.
firstdog = true; // 2.
// ||Dx sN||2
// Dx = diag((Sx)1, ... ,(Sx)n)
// se almacena como Sx E R^n..
// Para implementacion sugerida de expresiones
// que conciernen a Dx, ver seccion 1.5 del apendice
// 3.
Dx_sN = Dx.diagonalProduct(sN);
Newtlen = Dx_sN.norm2();
// Calcula y checa un nuevo paso
do { // 4.
// Encuentra el nuevo paso con el algoritmo double dogleg
// A.6.4.4
// respetando la declaracion
LWAIT();
LLOG << "\n------------------------------";
doubleDogLegStep(n, g, L,
sN, Sx, Newtlen,
maxstep,
delta, firstdog,
Cauchylen, eta, sSD,
v,
s, Newttaken); // 4.2
// Checa un nuevo punto y actualiza el radio
// de confianza A6.4.5.
// La segunda L en la lista de parametros es dummy
// test s[1] = -0.334;
// test s[2] = -0.335;
LLOG << "\nSe uso el paso de "
<< (Newttaken ? "Newton" : "Dogleg");
LWAIT();
LLOG << "\n------------------------------";
updateTrustRegion(n, xc, fc,
theFunction, g, L,
s, Sx, Newttaken,
maxstep, steptol, steptype,
H,
delta, retcode, xpprev,
fpprev,
xp, fp, maxtaken); // 4.3.
LLOG << "\ndelta = " << delta;
/*void updateTrustRegion(int n, VectorDouble &xc, double fc,
FunctionND_1D &theFunction, VectorDouble &g, SquaredMatrix &L,
VectorDouble &s, VectorDouble &Sx, bool Newttaken,
double maxstep, double steptol, int steptype,
SquaredMatrix &H,
double &delta, int &retcode, VectorDouble &xpprev,
double &fpprev,
VectorDouble &xp, double &fp, bool &maxtaken) {
*/
} while (retcode > 1); // 4U. until retcode < 2
// didac
LLOG << "\n-----IN/OUT DESPUES-----";
LLOG << "\ndelta c. Radio de la region de confianza. delta = "<< delta;
LLOG << "\n-----OUT-----";
LLOG << "\nBandera de retorno. 0: x+ satisf., 1: x+ no satisf. "
<< "\nretcode = " << retcode;
LLOG << "\nxp = " << xp;
LLOG << "\nfp = " << fp;
LLOG << "\nmaxtaken = " << maxtaken;
LWAIT();
LLOG << "\n } doubleDogLegDriver";
}
//----------------------------------------------------------------------------
#endif
/*
try {
// Presentacion
cout << "OptDogLeg. "
<< "Daniel Alba Cuellar, Omar Posada Villarreal";
// Entrada
cout << "\nSe usara logica 1.";
cout << "\nLa funcion vectorial es de dimension 2.";
// Var de entrada
int n = 2; //dimension
double fc; // escalar
VectorDouble xc(n, 1);
VectorDouble Sx(n, 1);
SquaredMatrix H(n, 1);
Function1 func;
/* Comentar al probar
cout << "\nPunto inicial Ejemplo: 1 1. xc: ";
xc.readElements();
cout << xc;
cout << "\nEscala Ejemplo: 1 1. Sx: ";
Sx.readElements();
cout << Sx;
/**/
/*
// Proceso
VectorDouble g(n, 1);
SquaredMatrix L(n, 1);
VectorDouble sN(n, 1);
//VectorDouble Sx(n, 1);
double maxstep;
double delta;
// Var doubleDogLegDriver
double steptol;
int retcode;
VectorDouble xp(n, 1);
double fp;
bool maxtaken;
/*
// Var solo de doubleDogLegStep()
double Newtlen;
bool firstdog;
double Cauchylen;
double eta;
VectorDouble sSD(n, 1);
VectorDouble v(n, 1);
VectorDouble s(n, 1);
bool Newttaken;
*/
/*
SquaredMatrix LT_L(n, 1);
SquaredMatrix LT_LInv(n, 1);
/* Comentar al entregar
/*
xc = [ 1 1 ]
H(xc) =
[ 14 0 ]
[ 0 2 ]
H^-1(xc) =
[ 0.0714286 0 ]
[ 0 0.5 ]
*/
/* xc[1] = 1;
xc[2] = 1;
Sx[1] = 1;
Sx[2] = 1;
H[1][1] = 14;
H[1][2] = 0;
H[2][1] = 0;
H[2][2] = 2;
if (!factorizeCholesky(n, H, L) ) {
throw runtime_error("No se pudo factorizar.");
}
g[1] = 6.0;
g[2] = 2.0;
//LT_L = L.transpose() * L;
LT_LInv[1][1] = 1.0/14.0;
LT_LInv[2][2] = 1.0/2.0;
//error sN = (LT_LInv * g) * (-1.0);
/* doglegstep*/
/* sN[1] = -0.428571;
sN[2] = -1.0;
/**/
/* driver
sN[1] = 0.334;
sN[2] = 0.335;
/**/
/* // Es una matriz cuadrada, pero se almacenara como un vector a Sx
VectorDouble Dx(Sx);
VectorDouble Dx_sN(n, 1);
Dx_sN = Dx.diagonalProduct(sN);
//Newtlen = Dx_sN.norm2();
// Dado por el usuario
maxstep = 1.0;
// Dado por el usuario
delta = 0.5; // dogleg 0.75; //ejemplo
//firstdog = true; // init updTrusReg
// Init solo doglegdriver
fc = func.evaluate(xc);
// Precision Dado por el usuario
steptol = 10E-6;
// Mostrar parametros
cout << "\nHessiana. H:\n" << H;
waitEnter();
cout << "\n-----IN-----";
cout << "\nDimension. n = " << n;
cout << "\nPunto inicial. xc: " << xc;
cout << "\nf(xc). Funcion en el punto inicial. fc: " << fc;
cout << "\ngrad(f(xc)). g = " << g;
cout << "\nTriangular inferior del Hessiano. L:\n" << L;
cout << "\nsc^N = -Hc^-1 grad(f(xc)). sN:" << sN;
cout << "\nEscala metrica. Sx: " << Sx;
cout << "\nMaxima delta. maxstep = " << maxstep;
cout << "\nsteptol = " << steptol;
cout << "\n-----IN/OUT ANTES-----";
cout << "\ndelta c. Radio de la region de confianza. delta = "<< delta;
waitEnter();
//cout << "\nPaso de Newton. Newtlen = " << Newtlen;
//cout << "\ntrue: 1a. iteracion. firstdog = " << firstdog;
cout << "\n------------------------------";
doubleDogLegDriver(n, xc, fc,
func, g, L,
sN, Sx,
maxstep, steptol,
delta,
retcode, xp, fp, maxtaken);
cout << "\n-----IN/OUT DESPUES-----";
cout << "\ndelta c. Radio de la region de confianza. delta = "<< delta;
//cout << "\n||s^CP||2. Cauchylen = " << Cauchylen;
//cout << "\nFactor para pasar de CP a ^N. eta = " << eta;
//cout << "\ns^CP, vector de xc a CP. sSD = " << sSD;
//cout << "\neta s^(^N) - s^CP. Vector para pasar de CP a ^N. v = " << v;
cout << "\n-----OUT-----";
cout << "\nBandera de retorno. 0: x+ satisf., 1: x+ no satisf. "
<< "\nretcode = " << retcode;
cout << "\nxp = " << xp;
cout << "\nfp = " << fp;
cout << "\nmaxtaken = " << maxtaken;
//cout << "\ns^CP + lambda v. Paso. s = " << s;
//cout << "\ntrue: se tomo el paso de Newton. Newttaken = "
// << toString(Newttaken);
} catch (exception e) {
cout << "\nmain: Excepcion " << e.what();
}
*/
// Fin------------------------------------------------------------------------