/** @(#) updatetrustregion.cpp */
// #include "ap.h"
/** @(#) .h */
/**
* @author Daniel Alba Cuellar
* @author Omar Posada Villarreal
* @version 1.0, /03/2002
*/
#ifndef UPDATETRUSTREGION_H
#define UPDATETRUSTREGION_H
//----------------------------------------------------------------------------
#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;
/** Dado un paso s, encontrado por el algoritmo Hook Step o DogStep,
* decide cuando x+ = xc + s debe ser aceptado para la siguiente iteracion.
* Si se elige, escoge la region de confianza inicial para
* la siguiente iteracion.
* Si no, disminuye o aumenta la region de confianza para la iteracion actual.
* 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 de la Hessiana
* @param s Paso: dogleg o hook
* @param Sx Escala. Valores usuales.
* @param Newttaken true: Se tomo el paso de Newton
* @param maxstep Maximo taman~o del paso, delta\. Dado por QuasiNewton
* @param steptol Precision del paso
* @param steptype
* 1: hookstep. Necesita a H
* 2: dogstep
* @param H Modelo Hessiano simetrica. Usada si steptype = 1.
* @param delta in/out Radio de la region de confianza.
* @param retcode in/out
* 0: x+ aceptado como la siguiente iteracion,
* delta es la region de confianza para la siguiente iteracion.
* 1: x+ no satisfactorio pero el paso global es terminado
* debido a que el taman~o relativo de (x+ - xc)
* esta dentro de la tolerancia del taman~o de paso para detener
* el algoritmo de detencion, indicando en la mayoria de los casos
* que no sera posible un progeso.
* Ver termcode = 3 de A7.2.1. y A7.2.3.
* 2: f(x+) es muy grande, la iteracion actual va a continuar con una
* delta reducida.
* 3: f(x+) suficientemente pequen~a, pero con la oportunidad de tomar
* un paso exitoso mas largo, que parece ser suficientemente
* bueno que el de la iteracion actual,
* el cual va a continuar con una nueva delta doble.
* @param xpprev in/out x+prev xpprev, fpprev tendran valores actuales en la
* entrada, o nyuevos valores en la salida solo si retcode = 3, en ese punto.
* Es la posicion siguiente aceptable antes de probar con una nueva x+.
* @param fpprev in/out f+prev. Es el valor de la funcion aceptable antes
* de probar con una nueva x+.
* @param xp return x+ Siguiente paso elegido.
* @param fp return f+ Valor de la funcion con el siguiente paso elegido.
* @param maxtaken return
* true: Si el taman~o de paso es mayor que el 99% de maxstep.
* Acepta x+ como una nueva iteracion, y escoge nueva delta.
* @version 1.0, 20/04/2002, sig vers LowerTriangularMatrix &L
* @see Dennis. A6.4.5. TRUSTREGUP pagina 338
*/
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,
int logType = SHORT_LOG, int waitType = NO_WAIT)
throw (invalid_argument) {
LLOG << "\n{ updateTrustRegion";
// Var locales de A6.4.3
VectorDouble Dx_s(n, 1);
// Es una matriz cuadrada, pero se almacenara como un vector de Sx
VectorDouble Dx(Sx);
double alfa;
double aux;
double delta_f; // incremento en f
double delta_fpred;
double deltatemp;
double initslope;
double mx; // maximo
double rellength; // Longitud relativa del paso
double steplen;
double sum;
double temp;
int i, j;
// Significa una llamada inicial a A6.4.4
// A6.4.5.
maxtaken = false; // 1.
// alfa es una constante usada en la prueba de aceptacion del paso
// f(x+) <= f(xc) + alfa g^T (x+ - xc).
//Puede ser cambiado al modificar este enunciado
alfa = 10E-4; // 2.
// 3.
// ||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
// ??? Dx = Sx.takeDiagonal();
Dx_s = Dx.diagonalProduct(s);
steplen = Dx_s.norm2();
xp = xc + s; // 4.
fp = theFunction.evaluate(xp); // 5.
delta_f = fp - fc; // 6.
initslope = g * s; // 7. innerProduct. g^T s
LLOG << "\nx+, xp = " << xp;
LLOG << "\nf(x+), fp = " << fp;
LLOG << "\ndelta_f = " << delta_f;
LLOG << "\ninitslope = " << initslope;
if (retcode != 3) { // 8.
// Puede ser necesario para prevenir un error runtime en la sig linea
fpprev = 0.0;
}
//9a.
if( (retcode == 3)
&&
( (fp >= fpprev)
|| (delta_f > (alfa * initslope))
)
) {
// Restaura x+ a x+prev y termina el paso global
retcode = 0; // 9a.1.
xp = xpprev; // 9a.2.
fp = fpprev; // 9a.3.
delta /= 2.0; // 9a.4.
LLOG << "\n --> updateTrustRegion 9a.";
// Esta comentado porque manejamos minimizacion,
// no buscamos las raices
//FVp = Fpprev; // 9a.5. Solo ecuaciones no lineales
// A6.4.5. de aqui regresa.
} else if (delta_f >= (alfa * initslope)) { // 9b.
// 9b.1.
// f(x+) es muy grande
//rellength = max(i=[1,n]){ |s[i]| / (max{ |x+[i]|, 1/Sx[i])}
mx = fabs(s[1]) / max( fabs(xp[1]), 1.0/Sx[1] );
for (i = 1; i <= n; ++i) {
aux = fabs(s[i]) / max( fabs(xp[i]), 1.0/Sx[i] );
if (aux > mx) {
mx = aux;
}
}
rellength = mx;
if (rellength < steptol) { // 9b.2.
// xp - xc muy chico, termina el paso global
retcode = 1; // 9b.2T.1.
xp = xc; // 9b.2T.2.
LLOG << "\n --> updateTrustRegion 9b.2";
// A6.4.5. de aqui regresa
} else { // 9b.2E
// Reduce delta, continua el paso global
retcode = 2; // 9b.2E.1
// 9b.2E.2
deltatemp = -(initslope * steplen)
/ (2.0 * (delta_f - initslope));
if (deltatemp < (0.1 * delta)) { // 9b.2E.3a.
delta *= 0.1;
} else if (deltatemp > (0.5 * delta)) { // 9b.2E.3b.
delta *= 0.5;
} else {
delta = deltatemp; // 9b.2E.3c.
}
LLOG << "\n --> updateTrustRegion 9b.2E.";
// A6.4.5 de aqui regresa
}
} else { // 9c.
// f(x+) suficientemente chico
// 9c.1-2 Asigna delta_fpred = g^T s + 1/2 s^T H s
delta_fpred = initslope; // 9c.1.
if (steptype == 1) { // 9c.2.
// hookstep, calcula s^T H s
LLOG << "\nHookstep sin implementar";
} else { // 9c.2E
// dogleg step, calcula s^T L L^T s
for (i = 1; i <= n; ++i) { // 9c.2E.1.
// 9c.2E.1.1.
sum = 0.0;
for (j = i; j <= n; ++j) {
// Checar con [][]
sum += L(j, i) * s(j);
}
temp = sum;
delta_fpred += square(temp) / 2.0; // 9c.2E.1.2
} // 9c.2E.1.
// 9c.3.
if( (retcode != 2)
&&
( ( fabs(delta_fpred - delta_f)
<= (0.1 * fabs(delta_f))
)
|| (delta_f <= initslope)
) && (Newttaken == false)
&& (delta <= 0.99 *maxstep)
) { // 9c.3T.
// doblar delta y continuar paso global
retcode = 3; // 9c.3T.1.
xpprev = xp; // 9c.3T.2.
fpprev = fp; // 9c.3T.3.
delta = min( (2.0*delta), maxstep); // 9c.3T.4.
LLOG << "\n --> updateTrustRegion 9c.3T.";
// Esta comentado porque manejamos minimizacion
// no buscamos las raices
// Solo ecuaciones no lineales
//Fpprev = FVp; // 9c.3T.5.
// A6.4.5 de aqui regresa
} else { // 9c.3E.
// Acepta x+ como una nueva iteracion,
// y escoge nueva delta
retcode = 0; // 9c.3E.1.
if (steplen > (0.99 * maxstep)) { // 9c.3E.2.
maxtaken = true;
}
// 9c.3E.3a.
if (delta_f >= (0.1 * delta_fpred)) {
// decrementa delta para sig iteracion
delta /= 2.0;
// 9c.3E.3b.
} else if (delta_f <= (0.75 * delta_fpred)) {
// incrementa delta para sig iteracion
delta = min( (2.0*delta), maxstep);
} // 9c.3E.3c else { no se cambia la delta }
} // 9c.3E.
} // 9c.2E.
} // 9c.
LLOG << "\ndelta_fpred = " << delta_fpred;
LLOG << "\n } updateTrustRegion";
} // A6.4.5.
//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------