downhill simplex

juri23

Hi habe bis jetzt nur einen Bergsteiger Algorithmus geschrieben. Nun wollte ich mich an einen Simplex rantasten. ( Habe einen fertigen gefunden )
simplex downhill Algorithmus :
Hier der Code : Würde mich sehr freuen wenn sich jemand meine Notizen anschaut ,ob ich mit meinen Vermutungen richtig liege. Ganz besonders in Bezug auf die 3 Punkte
P[1][1]= 1.0; P[1][2]=2.0;…………..
Außerdem wenn ich ein Maximum suche müsste ich doch nur die Vorzeichen ändern , also aus < wird >. Habe es auch so bei meinem selbstgeschriebenen Bergsteiger Alg. Gemacht. .lg

/**********************************************************
*   Multidimensional minimization of a function FUNC(X)   *
*  where X is an NDIM-dimensional vector, by the downhill *
*  simplex method of Nelder and Mead.                     *
* ------------------------------------------------------- *
* SAMPLE RUN: Find a minimum of function F(x,y):          *
*             F=Sin(R)/R, where R = Sqrt(x*x+y*y).        *
*                                                         *
*  Number of iterations: 22                               *
*                                                         *
*  Best NDIM+1 points:   ( Die Koordinaten  )             *
*   4.122686  1.786915                                    *
*   4.166477  1.682698                                    *
*   4.142454  1.741176                                    *
*                                                         *
*  Best NDIM+1 mimimum values:   //DER Y_WERT             *
*   -0.2172336265                                         *
*   -0.2172336281                                         *
*   -0.2172336271                                         *
*                                                         *
* ------------------------------------------------------- *
* REFERENCE: "Numerical Recipes, The Art of Scientific    *
*             Computing By W.H. Press, B.P. Flannery,     *
*             S.A. Teukolsky and W.T. Vetterling,         *
*             Cambridge University Press, 1986"           *
*             [BIBLI 08].                                 *
*                                                         *
*                      C++ Release By J-P Moreau, Paris.  *
**********************************************************/
#include <stdafx.h>
#include <iostream.h>
#include <stdio.h>
#include <math.h>

#define   MP  22   //Anzahl der Iterationen
#define   NP  21   //Maximum value for NDIM=20 //Anzahl der Parameter "Dimension des Parameterraums"

typedef   double MAT[MP][NP];

MAT P;
double Y[MP], PT[MP];
int I,ITER; //Anzahl der Iterationsschritte
int J,NDIM; // NDIM = Dimensionen
double FTOL; // Ist der Kleinste Funktionswert ? // Toleranzen

// user defined function to minimize
double FUNC(double *P) {
  double R;
  R=sqrt(P[1]*P[1]+P[2]*P[2]);
  if (fabs(R) < 1e-12)
    return 1.0;
  else
    return sin(R)/R;
}

void AMOEBA(MAT P, double *Y, int NDIM, double FTOL, int *ITER) {
/*-------------------------------------------------------------------
! Multidimensional minimization of the function FUNC(X) where X is
! an NDIM-dimensional vector, by the downhill simplex method of
! Nelder and Mead. Input is a matrix P whose NDIM+1 rows(Zeilen) are NDIM-
! dimensional vectors which are the vertices ( Ecken )of the starting simplex
! (Logical dimensions of P are P(NDIM+1,NDIM); physical dimensions
! are input as P(NP,NP)). Also input is the vector Y of length NDIM
! +1, whose components must be pre-initialized to the values of FUNC
! evaluated at the NDIM+1 vertices (rows) of P; and FTOL the fractio-
! nal convergence tolerance to be achieved in the function value. On
! output, P and Y will have been reset to NDIM+1 new points all within
! FTOL of a minimum function value, and ITER gives the number of ite-
! rations taken.

Zusammengefasst:
- In der Funktion FUNC(x) soll das Minimum gefunden werden.
- X ist ein N-Dimensioaler Vektor.
- Y ist der Eingangsverktor
- Als eingangswert dient die Matrix P mit Dimension+1 Zeilen welche die Eckt Punkte ( Anfangswerte) beschreibt.
- P(NP,NP) beschreibt die Anzahl der Parameter

!-------------------------------------------------------------------*/
// Label:  e1
const int NMAX=20, ITMAX=500; //Maximale Anzahl von Dimensionen und Interationen

//Expected maximum number of dimensions, three parameters which define
// the expansions and contractions, and maximum allowed number of
//iterations.
  double PR[MP], PRR[MP], PBAR[MP];
  double ALPHA=1.0, BETA=0.5, GAMMA=2.0;
  int I,IHI,ILO,INHI,J,MPTS;
  double RTOL,YPR,YPRR;
  MPTS=NDIM+1;
  *ITER=0;
e1:ILO=1;
  if (Y[1] > Y[2]) {
    IHI=1;
    INHI=2;
  }
  else {
    IHI=2;
    INHI=1;
  }
  for (I=1; I<=MPTS; I++) {
    if (Y[I] < Y[ILO])  ILO=I;
    if (Y[I] > Y[IHI]) {
      INHI=IHI;
      IHI=I;
    }
    else if (Y[I] > Y[INHI])
      if (I != IHI)  INHI=I;
  }
//Compute the fractional range from highest to lowest and return if
//satisfactory.
  RTOL=2.0*fabs(Y[IHI]-Y[ILO])/(fabs(Y[IHI])+fabs(Y[ILO]));
  if (RTOL < FTOL)  return;  //normal exit
  if (*ITER == ITMAX) {
    printf(" Amoeba exceeding maximum iterations.\n");
    return;
  }
  *ITER= (*ITER) + 1;
  for (J=1; J<=NDIM; J++)  PBAR[J]=0.0;
  for (I=1; I<=MPTS; I++)
    if (I != IHI)
      for (J=1; J<=NDIM; J++)
        PBAR[J] += P[I][J];
  for (J=1; J<=NDIM; J++) {
    PBAR[J] /= 1.0*NDIM;
    PR[J]=(1.0+ALPHA)*PBAR[J] - ALPHA*P[IHI][J];
  }
  YPR=FUNC(PR);
  if (YPR <= Y[ILO]) {
    for (J=1; J<=NDIM; J++)
      PRR[J]=GAMMA*PR[J] + (1.0-GAMMA)*PBAR[J];
    YPRR=FUNC(PRR);
    if (YPRR < Y[ILO]) {
      for (J=1; J<=NDIM; J++) P[IHI][J]=PRR[J];
      Y[IHI]=YPRR;
    }
    else {
      for (J=1; J<=NDIM; J++) P[IHI][J]=PR[J];
      Y[IHI]=YPR;
	}
  }
  else if (YPR >= Y[INHI]) {
	if (YPR < Y[IHI]) {
      for (J=1; J<=NDIM; J++)  P[IHI][J]=PR[J];
      Y[IHI]=YPR;
    }
    for (J=1; J<=NDIM; J++)  PRR[J]=BETA*P[IHI][J] + (1.0-BETA)*PBAR[J];
    YPRR=FUNC(PRR);
    if (YPRR < Y[IHI]) {
      for (J=1; J<=NDIM; J++)  P[IHI][J]=PRR[J];
      Y[IHI]=YPRR;
    }
    else
      for (I=1; I<=MPTS; I++)
		if (I != ILO) {
		  for (J=1; J<=NDIM; J++) {
            PR[J]=0.5*(P[I][J] + P[ILO][J]);
	        P[I][J]=PR[J];
		  }
          Y[I]=FUNC(PR);
		}
  }
  else {
    for (J=1; J<=NDIM; J++)  P[IHI][J]=PR[J];
    Y[IHI]=YPR;
  }
  goto e1;
}

void main()  {

  NDIM=2;      // Anzahl der Variablen
  FTOL=1e-8;   // tolerance

  //define NDIM+1 initial vertices (one by row) // Die Eckpunkte werden hier deferiert ( Alles Vermutungen )
  P[1][1]= 1.0; P[1][2]=2.0;    //Pro Zeile immer nur einer wir haben 2 Variablen (2 Dimensionen) also müssen wir 3 Punkte Definieren
  P[2][1]=-2.0; P[2][2]=-3.0;   //Wenn ich richtig liege, gibt man für jeden Parameter zwei verschiedene Werte vor, berechnet den ersten Punkt
  P[3][1]= 4.0; P[3][2]=2.0;    // mit allen Erstwerten und die anderen NPunkte auch mit den Erstwerten, nur das jeweils einer der Parameter durch seinen 
                                // zweitwert ersetzt wird.

  // Wir berechnen an jedem Punkt 1,2,-2,-3,4,2 die Funktion und dadurch wird unser Kordinatensatz verändert
  // Der eine Wert beschreibt die Koordinate der andere seinen Funktionswert
  // Oder ist er nur ein Funktionswert erweitert um eine Dimension 

//Initialize Y to the values of FUNC evaluated // In die Y Martix packen wir unsere Anfangswerte
//at the NDIM+1 vertices (rows] of P
  for (I=1; I<=NDIM+1; I++) {
    PT[1]=P[I][1]; PT[2]=P[I][2];
    Y[I]=FUNC(PT);
  }

  //call main function
  AMOEBA(P,Y,NDIM,FTOL,&ITER);

  //print results
  printf("\n Number of iterations: %d\n\n", ITER);
  printf(" Best NDIM+1 points:\n");
  for (I=1 ; I<=NDIM+1; I++) {
    for (J=1; J<=NDIM; J++)  printf("  %f", P[I][J]);
    printf("\n");
  }
  printf("\n Best NDIM+1 mimimum values:\n");
  for (I=1; I<=NDIM+1; I++) printf(" %14.10f\n", Y[I]);
  printf("\n");

  system( "pause" );

}

// end of file tamoeba.cpp

juri23

ich hoffe hier im Forum ist jemand der sich die mühe kurz macht. Würde mich sehr freuen

Tu dir einen Gefallen und such eine andere Implementierung. Diese hier ist einfach nur grauenhaft.

juri23

meinst du damit meine Notizen ?
//...

das wird so nicht bleiben sind wie Notizzettel.

Nein, nicht die Notizen. Ich meine den Code den du gefunden hast. Dich trifft keine Schuld

juri23

du meinst wegen goto ?
nobody is perfect.
Aber trotzdem könntest du mir weiter helfen ?

otze

worauf wir raus wollen:

double PR[MP], PRR[MP], PBAR[MP];
  double ALPHA=1.0, BETA=0.5, GAMMA=2.0;
  int I,IHI,ILO,INHI,J,MPTS;
  double RTOL,YPR,YPRR;

Das dechiffriert hier niemand.

juri23

http://de.wikipedia.org/wiki/Downhill-Simplex-Verfahren
wiki hat dazu was interessantes.

// Deklarationen

// Benutzervorgaben, -eingaben
int NP= ... ; // Anzahl der Parameter, die in FQS() eingehen,
// auch "Dimension des Parameterraums"

double ParP[NP]={ ... }; // primärer Parametersatz
double ParS[NP]={ ... }; // sekundärer Parametersatz, gibt Schwankungsbreite vor

double FQS(int pf) // Zu minimierende Funktion, hängt
{ // von den Parametern Par0[pf][...] ab.
... // Muss je nach Problemstellung erstellt werden.
// So schreiben, dass sie für den Parametersatz Nr. pf
} // berechnet und am Schluss in Par0[pf][0] (s. u.)
// abgespeichert wird.

// Ende Benutzervorgaben

// Initialisierung

// Simplex-Parameter
double NA = 1, NB = .5, NC = 2; // Alpha, Beta, Gamma

// Parametersätze, NP+1 Stück für den Simplex selbst, plus 3 weitere
double Par0[NP+4][NP+1];
// darin 1. Index: 0..NP die (NP+1) Parametersätze des Simplex,
// dabei auf 0 immer der schlechteste Punkt des Simplex,
// auf 1 immer der beste Punkt des Simplex
// NP+1 Punkt P*, s. u.
// NP+2 Punkt P**, s. u.
// NP+3 Mittelpunkt P-Quer des Simplex, s. u.
// darin 2. Index: 0 Fehlerwert für diesen Satz
// 1..NP einzelne Parameterwerte dieses Satzes

// Initialisierung

// Simplex aus Primär- und Sekundärwerten
Parametersätze 0 bis NP mit Primärwerten belegen;
In Parametersätzen 1 bis NP (i) jeweils
Parameter Nr. i mit Sekundärwert belegen;

// Berechnung der Anfangs-Fehlersummen
für Parametersätze 0 bis NP: feh=FQS(i); // speichert Ergebnis in Par0[i][0], s. o.

// Schleife für Iteration, Hauptschleife
bool fertig = false;
while (!fertig)
NT++; // Simplex-Schritt Nr. NT

// schlechtesten und besten Punkt ermitteln (höchster bzw. niedrigster Fehlerwert)
Werte fmin und fmax auf Fehlerwert von Satz Nr. 0 vorbelegen;
für alle Parametersätze von 1 bis NP:
wenn Fehler für Satz Nr. i größer als fmax:
noch schlechterer Punkt, fmax auf diesen Wert, Punkt mit Nr. 0 tauschen
wenn Fehler für Satz Nr. i kleiner als fmin:
noch besserer Punkt, fmin auf diesen Wert, Punkt mit Nr. 1 tauschen

// Konvergenzabfrage, sind wir fertig?
// Dies muss auch je nach Problem individuell erstellt werden,
// hier als Beispiel die Bedingung, dass alle Parameter nur noch um 1 Promille
// schwanken.

fertig = true;
für alle Parameter in allen Parametersätzen von 0 bis NP:
wenn nur für einen Parameter die relative Abweichung des Wertes
zwischen schlechtestem und bestem Punkt größer als 0.001, dann:
fertig=false; // doch noch nicht fertig

//Wenn fertig, dann beende mit bestem Punkt Par0[1][1..3] als Ergebnis des Algorithmus

// Mittelpunkt P-QUER des Simplex nach Index NP+3
Mittelwerte über Sätze 1 bis NP (also ohne schlechtesten Punkt Nr. 0!) bilden
und als Satz Nr. NP+3 speichern;

// Punkt P* durch Reflexion des schlechtesten Punkts
// am Mittelpunkt P-QUER (mit Nelder/Mead-Parameter alpha=NA) nach Index NP+1
innerhalb eines Parametersatzes:
Par0[NP+1][i]=(1 + NA) * Par0[NP+3][i] - NA * Par0[0][i];
fs=Par0[NP+1][0]=FQS(NP+1); // und Fehler berechnen

// Fallunterscheidung, ob Verbesserung erreicht
wenn fs<fmin // neuer Punkt P* im Vergleich zu bisher bestem
/ // ja, besser!
| // weiteren Schritt zu Punkt P** in dieselbe Richtung nach Index NP+2
| // (mit Nelder/Mead-Parameter gamma=NC)
| innerhalb eines Parametersatzes:
| Par0[NP+2][i]=(1 + NC) * Par0[NP+1][i] - NC * Par0[NP+3][i];
| fs=Par0[NP+2][0]=FQS(NP+2); // und Fehler berechnen
|
| wenn fs>=fmin // wieder Vergleich mit bisher bestem Punkt
| / // keine weitere Verbesserung: P** verwerfen, P* nach Index 0,
| | // also bisher schlechtesten Wert durch neuen (etwas besseren) ersetzen
| \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren;
| else // von fs>=fmin
| / // ja, auch eine Verbesserung!
| | // auch besser als P* ?
| | wenn fs>Par0[NP+1][0]
| | / // nein, P* weiterverwenden und damit bisher schlechtesten auf Index 0 ersetzen
| | \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren;
| | else // von fs>Par0[NP+1][0]
| | / // ja, P** weiterverwenden und damit bisher schlechtesten auf Index 0 ersetzen
\ \ \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren;

else // von fs<fmin
/ // nicht besser als bisher bester, nun prüfen, ob P* wenigstens überhaupt eine Verbesserung
| wenn fs von P* für wenigstens einen der Punkte 0..NP kleiner ist, dann:
| / // ja, besser als mindestens einer, mit P* den bisher schlechtesten Punkt auf Index 0 ersetzen
| \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren;
| else // von fs von P*
| / // nein, weiterhin schlechter als die anderen Punkte
| | // Zusatzprüfung: sogar schlechter als bisher schlechtester Punkt fmax?
| | wenn fs>fmax // schlechter als bisher schlechtester?
| | // ja, erstmal nichts tun
| | else // von fs>fmax
| | / // nein, also wenigstens bisher schlechtesten Punkt damit ersetzen
| \ \ Parametersatz Nr. NP+1 nach Nr. 0 kopieren;
| // neuer Punkt P** durch Kontraktion zwischen bisher schlechtestem Punkt (Index 0)
| // und P-QUER (Index NP+3) nach Index NP+2 (mit Nelder/Mead-Parameter beta=NB)
| innerhalb eines Parametersatzes:
| Par0[NP+2][i]= NB * Par0[0][i] + (1 - NB) * Par0[NP+3][i];
| fs=Par0[NP+2][0]=FQS(NP+2); // und Fehler berechnen
| // P** besser als bisher schlechtester Punkt?
| wenn fs<fmax // besser als bisher schlechtester?
| / // ja, bisher schlechtesten durch P** ersetzen
| \ Parametersatz Nr. NP+2 nach Nr. 0 kopieren;
| else // von fs<fmax
| / // nein, keinerlei Verbesserung, Notfall
| | // folgt allgemeine Kompression um Faktor 2 um den bisher besten Punkt herum
| | für alle Parametersätze außer Nr. 1 (bisher bester): // also j=0 und 2..NP
| | innerhalb eines Parametersatzes j:
| | Par0[j][i]=(Par0[j][i] + Par0[1][i]) / 2;
\ \ fs=Par0[j][0]=FQS(j); // und Fehler berechnen
// Schleifenende