/** @(#) numericalap/logisticnewtondogleg.cpp */
#include "logisticnewtondogleg.h"
class LeastSquareFunction : public FunctionND_1D {
private:
int m;
VectorDouble t;
VectorDouble P;
public:
LeastSquareFunction(VectorDouble &time, VectorDouble &population)
throw (invalid_argument) {
if (time.getSize() == population.getSize()) {
m = time.getSize();
} else {
throw invalid_argument(
"LeastSquareFunction: Dimensiones distintas");
}
t = time;
P = population;
}
// No se checa las dimensiones de xc y Fc.
double evaluate(VectorDouble &xc) {
xc.setLogic(1);
// Funcion del usuario
// f(x1, x2) = x1^4 + x1^2 + x2^2
// xc = [1 , 1]
double a = xc[1];
double b = xc[2];
double c = xc[3];
double fc = 0.0; // ==sum
int i;
// Residuals
// manualmente 0.5 (R^t R)
for (i = 1; i <= m; ++i) {
fc += square(
( a * c / ( b * c + exp(-a * t[i]) )
) - P[i]
);
}
fc /= 2.0;
return fc;
}
}; // No olvidar ";"
/**
* @param t Tiempo. Dim: m.
* @param P Poblacion. Dim: m.
* @param x Coeficientes. x[1] = a, x[2] = b, x[3] = c. Dim: 3.
* @param ri return Residuos. Dim: m.
* @param g return Gradiente. Dim: m.
* @param J return Jacobiana. Dim: mx3.
* @param H return Hessiana. Dim: mxm.
*/
void logisticMatrices(VectorDouble &t, VectorDouble &P, VectorDouble &xc,
VectorDouble &ri, RectangularMatrix &J,
VectorDouble &g, SquaredMatrix &H) {
// Coeficientes
double a = xc[1]; // a: crecimiento (maltusiano),
double b = xc[2]; // b: competencia,
double c = xc[3]; // c: parametro de integracion
//(1er punto de la curva)
int m = t.getSize();
int i;
// Informacion de las 2as derivadas de la expansion de Taylor
SquaredMatrix S(3, 1, true); // inicializar a ceros
// Residuals
for (i = 1; i <= m; ++i) {
ri[i] = ( a * c / (b * c + exp(-a * t[i])) ) - P[i];
}
// Jacobiano
double eati, den2, den3; // aux
for (i = 1; i <= m; ++i) {
eati = exp(-a * t[i]);
den2 = square((b * c) + eati);
J(i, 1) = ( (b * square(c)) + ((1 + a*t[i]) * c * eati) ) / den2;
J(i, 2) = (-a * square(c)) / den2;
J(i, 3) = (a * eati) / den2;
}
// Gradiente = J^t ri
RectangularMatrix Jt(3, m, 1, false);
Jt.transpose(J);
g = Jt * ri;
// Logica cero
// Hessiana
// grad^2 ri. m matrices 3x3
SquaredMatrix *pHR = new SquaredMatrix[m];
for (i = 0; i < m; ++i) {
//ERROR
pHR[i].resize(3, 1, false);
}
for (i = 0; i < m; ++i) {
eati = exp(-a * t[i+1]);
den3 = pow((b * c) + eati, 3.0);
pHR[i](1, 1) =
( (-c * t[i+1] * eati) *
( (a * b * c * t[i+1]) - (a * t[i+1] * eati)
- (2 * b * c) - (2 * eati)
)
) / den3;
// factorizar ( -square(c) * (b * c + eati + 2 * a * t[i] * eati))
pHR[i](1, 2) =
( -square(c) * (b * c + eati + 2 * a * t[i+1] * eati)
) / den3;
pHR[i](1, 3) =
( -eati *
( (a * b * c * t[i+1])
- (a * t[i+1] * eati)
- (b * c) - eati
)
) / den3;
pHR[i](2, 2) = (2 * a * c*c*c) / den3;
pHR[i](2, 3) = (-2 * a * c * eati) / den3;
pHR[i](3, 3) = (-2 * a *b * eati) / den3;
}
// Calcular S. H = J^T J + S
SquaredMatrix sum(3, 1, true); // sumando = 0
for (i = 0; i < m; ++i) {
// sum = pHR[i] * ri[i];
SquaredMatrix::multiply(pHR[i], ri[i], sum);
//S = S + sum;
SquaredMatrix::add(S, sum, S);
}
// De triangular a simetrica
S.upperTriangularToSymmetric();
// H = J^T J + S
//H = J.transpose() * J + S;
RectangularMatrix aux(1, 1, false);
RectangularMatrix aux2(1, 1, false);
aux2.transpose(J);
RectangularMatrix::multiply(aux2, J, aux);
RectangularMatrix::add(aux, S, H);
}
/** 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 logisticNewtonDogLeg(int n, VectorDouble &t, VectorDouble &P,
VectorDouble &xc, double &fc,
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{ logisticNewton";
// log file
ofstream outFile(pathFile.data(), ios::out);
// Recordar cerrar outFile.close();
if ( !outFile ) { // !sobrecargado
throw runtime_error(
"logisticNewtonDogLeg: Error de apertura de archivo.");
}
// { posada
xc.setLogic(1);
Sx.setLogic(1);
//H.setLogic(1);
// } posada
SquaredMatrix H(n, 1); // Hessiana
SquaredMatrix QH(n, 1); // Q de 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))
// { logistic
int m = t.getSize();
RectangularMatrix J(m, 3, 1); // Jacobiano logistic
VectorDouble ri(m, 1); // Residuos
LeastSquareFunction theFunction(t, P);
SquaredMatrix eigenVectors(n, 1, true); //mu
VectorDouble initial(n, 1); // Para iniciar el eigen de la potencia
VectorDouble eVector(n, 1); // Eigen vector
initial.fill(1.0);
eVector.fill(1.0);
// } logistic
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);
double mu; // mu = -minEigVal(H)
SquaredMatrix mu_I(3, 1, false);
// newton iteration
for (k=1; k<=iterations; k++) {
//ERROR EN LA 2A CORRIDA
logisticMatrices(t, P, xc,
ri, J,
g, H);
// qrdcmp cambia a H, se toma una copia de H
QH.copyFrom(H);
// QH no es simetrica
sing = qrdcmp(QH, n, c, d); // aproximacion
// Si es singular, no se puede buscar la condicion.
if (sing) {
throw runtime_error(
"\nlogisticNewtonDogLeg: La matriz es singular.");
} else {
//sN Solucion aproximada a H(f(xk)) sk^N = -grad(f(xk))
R.createRMatrix(QH, d);
cond = R.condition();
if (cond > (1.0 / eta)) {
throw runtime_error(
"\nfiniteNewtonDogLeg: Hessiano mal condicionado.");
}
sN = g * -1.0;
qrsolv(QH, n, c, d, sN); //sN in/out
// Obtener triang inf de hessiana
if ( !factorizeCholesky(n, H, L)) {
// Perturbar la matriz.
//H = H + mu I, mu > 0
// positivo * 101%
// asegurarnos que no sea igual al min eigVal
/* Por Gerschgori*/
mu = -H.findMinLimitEigenValue() * (1.01);
//1 .05 = 1+ SQRT_MACHEPS);
/*
mu = fabs(H.dominantEigenValue(
initial, iterations, 1, eigenVectors, eVector)
) * (1.01);
/**/
// LA MU NO FUNCIONA
mu_I.toIdentity();
RectangularMatrix::multiply(mu_I, mu, mu_I);
RectangularMatrix::add(H, mu_I, H);
if ( !factorizeCholesky(n, H, L) ) {
throw runtime_error(
"logisticNewtonDogLeg: No se factorizo Hessiana.");
}
}
// 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
//antes
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.
*/
SLOG << "\t" << fc << "\t" << n << "\t" << xc;
FLOG << "\t" << fc << "\t" << n << "\t" << xc;
//if (retcode == 0) { // No hacer nada
if (retcode == 1) {
LLOG << "\n } logisticNewtonDogLeg";
outFile.close();
return;
}
// 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(
"\nlogisticNewtonDogLeg: 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 } logisticNewtonDogLeg";
outFile.close();
return;
}
SWAIT();
} // if
}
LLOG << "\n } logisticNewtonDogLeg";
outFile.close();
}
/** Checa parametros: dimensiones, rangos. */
void checkLogisticNewtonDogLeg(int n, VectorDouble &t, VectorDouble &P,
VectorDouble &xc, double &fc,
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("\nlogisticNewtonDogLeg: arg invalidos.");
}
// Prevee division entre cero
if (Sx.hasAlmostZeros()) {
throw invalid_argument(
"logisticNewtonDogLeg: Division entre cero.");
}
logisticNewtonDogLeg(n, t, P,
xc, fc,
Sx,
iterations,
maxstep, steptol,
delta,
eta,
logType, waitType,
pathFile);
}
// Fin-------------------------------------------------------------------------