brentmax.cc

//L'algoritmo restituisce il primo max relativo che trova calcolandolo con un errore relativo tol, guardando prima a destra dell'intervallo iniziale e poi a sinistra.

#include <iostream>
#include <cmath>

double goldmax(double func(double), double ax, double bx, double cx, double tol);
double expsearch(double func(double), double a, double b, double c, double tol, int k, double a0, double b0);

//Brentmax effettua una prima sezione aurea di [a,c] e passa i tre punti trovati a goldmax (se sono già 'buoni') oppure a expsearch.

//>>------------------------------
double brentmax(double func(double), double a, double c, double tol){
float b=c-((3.0-sqrt(5))/2)*(c-a);
if ((func(a) < func(b) && func(b) < func(c)) || (func(a) > func(b) && func(b) > func(c))){
return expsearch(func, a, b, c, tol, 0, a, b);
}else{
return goldmax(func, a, b, c, tol);
}
}
//>>-------------------------------


//Goldmax applica il metodo della sezione aurea per trovare il massimo a partire da un intervallo a<b<c tale che f(a)>f(b), f(c)>f(b). (questo codice è più o meno tradotto da quello che ho trovato all'indirizzo http://linneus20.ethz.ch:8080/1_5_2.html, da dove ho preso anche l'idea per il metodo esponenziale di brent)
double goldmax(double func(double), double ax, double bx, double cx, double tol){
double C=(3.0-sqrt(5))/2.0;
double R=1-C;

double x0=ax;
double x1;
double x2;
double x3=cx;

if (fabs(cx-bx) > fabs(bx-ax)){
x1=bx;
x2=bx+C*(cx-bx);
}else{
x2=bx;
x1=bx-C*(bx-ax);
}
double f1=func(x1);
double f2=func(x2);

int k=0;
while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))){ //sul modo in cui questa condizione tiene d'occhio l'errore relativo ho qualche vaga idea, ma non sono troppo sicuro.. però funziona..

    //std::cout << "k=" << k << "|a-b|=" << fabs(x3-x0) <<"<-->";
    if (f2 > f1){
x0=x1;
x1=x2;
x2=R*x1+C*x3; // x2=x1+C*(x3-x1)

      f1=f2;
f2=func(x2);
}else{
x3=x2;
x2=x1;
x1=R*x2+C*x0; // x1=x2+C*(x0-x2)
      f2=f1;
f1=func(x1);
}
++k;
}

double xmax;
if(f1 > f2){
xmax=x1;
}else{
xmax=x2;
}
return xmax;
}





//Expsearch è una funzione ricorsiva. Prepara a goldmax un intervallo adatto, partendo da un intervallo qualunque. Se quello fornito da fuori non va bene, comincia dilatandolo verso destra (metodo della ricerca esponenziale di brent); se dopo 100 passaggi non ha trovato un massimo in questa direzione, torna all'intervallo di partenza e ripete la procedura verso sinistra. Se non trova niente nemmeno qui fornisce 1234567890 come valore 'di errore'. Non sono riuscito a farle dare Inf come output..
double expsearch(double func(double), double a, double b, double c, double tol, int k, double a0, double b0){
if ((func(a) <= func(b) && func(b) <= func(c)) || (func(a) >= func(b) && func(b) >= func(c))){
double R=2/(sqrt(5)-1);
++k;
if(k<100){
//std::cout << "Expsearch (destra): livello " << k << "<-->";

      double x=c+R*(c-b);
return expsearch(func, b, c, x, tol, k, a0, b0);
}else if(k==100){                       
//std::cout << "Expsearch (sinistra, primo): livello " << k << "Valore" << "<-->";

      double x=a0-R*(b0-a0);
return expsearch(func, x, a, b, tol, k, a0, b0);
}else if(k<200){
//std::cout << "Expsearch (sinistra): livello " << k << "Valore" << a << "<-->";

      double x=a-R*(b-a);
return expsearch(func, x, a, b, tol, k, a0, b0);
}else{
return 1234567890;
}
}
return goldmax(func,a,b,c,tol);
}