function [lsg,iteration,nfc]=neldermead(fktname,x,iterationmax,c,subroutine)
% ======================================================================
% Eine Implementation der Polyeder-Methode von Nelder-Mead.
%
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% Helge Baumann, 11.8.1998
% ----------------------------------------------------------------------
% Aurufparameter:
%  (lediglich die ersten zwei Parameter muessen zwingend angegeben werden.)
% fktname : Name der Funktion F, deren Minimum bestimmt 
%   werden soll. An <fktname> wird ein Spaltenwektor(!) uebergeben. 
%   Es wird eine reelle Zahl als Bild erwartet.
% x : Spaltenvektor der Dimension n oder eine n x (n+1) Matrix.
%   Bei Vorgabe nur eines Spaltenvektors wird ein regelmaessiges 
%   Startpolyeder erzeugt. In diesem Fall muss dem Parameter c 
%   ein Wert uebergeben werden.
% iterationmax : Anzahl der maximalen Iterationen. Standard: 100
% c : Bei Vorgabe nur eines Spaltenvektors wird ein regelmaessiges
%   Startpolyeder gebildet. Die Variable c gibt dabei den Abstand je
%   zweier Startpunkte an. Standard: 1
% subroutine : Name einer Unterroutine, die in jedem Iterationsschritt
%   aufgerufen wird. Dient zum Protokolieren des Iterationsprozesses.
%   Wird <''> uebergeben, erfolgt kein Aufruf. Standard: ''
% ---------------------------------------------------------------------
% Rueckgabeparameter:
% lsg : Loesung (Spaltenvektor!)
% iteration : Anzahl der benoetigten Iterationen
% nfc : Anzahl der Funktionsauswertungen, muss allerdings in der
%   Unterroutine <fktname> gesetzt werden!
% ---------------------------------------------------------------------
% Programmvariablen:
% n  : Dimension der Vektorraumes
% m  : Anzahl der uebergebenen Vektoren (entweder m=1 oder m=n+1)
% x  : n x (n+1) Matrix, die n+1 Spaltenvektoren bestimmen das aktuelle
%      Polyeder
% fx : Vektor der Funktionswerte der n+1 Polyedervektoren
% b  : Index des Vektors mit dem "b"esten Funktionswert
% s  : Index des Vektors mit dem "s"chlechtesten Funktionswert
% z  : Index des Vektors mit dem "z"weitschlechtesten Funktionswert
% alpha, beta, gamma : Parameter der Polyeder-Methode
% epsilon : Abbruchgenauigkeit
% iterationmax : Anzahl der maximalen Iterationen
% abbruchbed,p2,p3,i,tt,tt2
%    :  Hilfsvariable
% ======================================================================
global nfc;

% Initialisierung
epsilon = 1E-10; alpha = 1; beta = 2; gamma = 0.5;
abbruchbed = 2*eps + 1;
iteration = 0;
nfc = 0;
lsg = [];

% --------------------------------------------------------------
% Hier findet eine Ueberpruefung der uebergebenen Parameter statt.
% Der eigentliche Algorithmus beginnt weiter unten.
error(nargchk(2,5,nargin))
if nargin < 5, subroutine = ''; end
if nargin < 4, c = 1; end
if nargin < 3, iterationmax = 10000; end

% Groesse der uebergebenen Vektoren bestimmen:
[n,m] = size(x);
% Ueberpruefen, ob die vorgegebenen Vektoren, die richtige Dimension haben.
if (n == 1) & (m > 2)
   disp('Bitte geben Sie einen Spaltenvektor vor!');
   return;
elseif (n == m+1) & (m > 1)
   disp('Bitte geben Sie Spaltenvektoren vor!');
   return;
elseif ~(n+1 == m) & (m >= 2)
   disp('Bitte geben Sie genau n+1 Startvektoren vor');
   return;
elseif (n > 1) & (m == 1)
   % wenn nur ein Vektor vorgegeben wurde, werden die anderen bestimmt:
   
   % Speicherplatz reservieren:
   x = cat(2,x,eye(n));
   
   % Parameter bestimmen:
   p2 = c * (sqrt(n+1)-1) / (sqrt(2)*n);
   p3 =  c /sqrt(2);
   
   % Startvektoren aufbauen:
   for i=2:n+1
      x(1:n,i) = x(1:n,1) + p2;
      x(i-1,i) = x(i-1,i) + p3;
   end;
end;
% Ende der Parameterueberpruefung
% --------------------------------------------------------------
% Beginn des Algorithmus
% --------------------------------------------------------------

% In x sind nun n+1 Spaltenvektoren der Dimension n gespeichert.
% Berechnung der Funktionswerte aller Vektoren:
[fx] = feval(fktname,x); 

% Subroutine aufrufen zu Protokollzwecken:
if ~isempty(subroutine), feval(subroutine,x,n); end


% Beginn der Iterationsschleife
% Beim ersten Durchgang ist abbruchbed immer groesser als epsilon
while (iteration < iterationmax) & (abbruchbed > epsilon)
   % Bestimmung des Index des besten Punktes
   [tt,b]=min(transpose(fx));
   % Bestimmung des Index des schlechtesten Punktes
   [tt,s]=max(transpose(fx));
   % Bestimmung des Index des zweitschlechtesten Punktes
   [tt,z] =max(fx([1:s-1,s+1:n+1]));
   if z>= s, z=z+1; end
   
   % Bestimme das Zentrum der n besseren Punkte:
   xquer = (sum(x(:,[1:s-1,s+1:n+1])')' ) / n;
   
  
   % Reflektiere den schlechtesten Punkts an xquer:
   xr = xquer + alpha * (xquer - x(:,s));
   [fxr] = feval(fktname,xr); 

   
   % Unterscheidung von 3 Faellen:
   
   % Ist f(xr) zwischen bestem und zweitschlechtestem Wert, so wird 
   % der schlechtest Punkt xs durch xr ersetzt.
   if (fx(b) <= fxr) & (fxr <= fx(z))
      x(:,s) = xr; fx(s) = fxr;
   % Ist f(xr) sogar noch besser als der bislang bekannte beste Wert...
   elseif (fxr < fx(b))
      % ...so koennte in der gleichen Richtung noch mehr zu holen sein.
      xe = xquer + beta*( xr - xquer);
      [fxe] = feval(fktname,xe); 
      % Ist die weitere Expansion erfolgreich, so tausche xe ansonsten xr
      % gegen den schlechtesten Punkt xs aus.
      if (fxe < fxr)
         x(:,s) = xe; fx(s) = fxe;
      else
         x(:,s) = xr; fx(s) = fxr;
      end;
      b = s; % auf dem Index des schlechtesten ist nun der beste Wert abgelegt.
   % Ist f(xr) schlechter als der zweitschlechteste Wert ist vermutlich
   % das Polyeder zu gross und eine Kontraktion erforderlich.
   elseif (fx(z) < fxr)
      if fxr >= fx(s)
         xk = xquer + gamma*( x(:,s)-xquer);
      else
         xk = xquer + gamma*( xr - xquer);
      end;
      [fxk] = feval(fktname,xk); 
      % Ist f(xk) besser als der schlechteste Wert, so wird xs gegen xk
      % ausgetauscht. Andernfalls wird das Polyeder um den besten Punkt
      % zusammengezogen.
      if (fxk < fx(s))
         x(:,s) = xk; fx(s) = fxk;
      else
         for i=1:n
            x(i,[1:b-1,b+1:n+1]) = (1 - gamma) * x(i,[1:b-1,b+1:n+1]) + gamma * x(i,b);
         end
         % Neuberechnung der Funktionswerte aller Vektoren nach Kontraktion:
         [fx] = feval(fktname,x); 
         
      end;
   end;
   
   % Subroutine aufrufen zu Protokollzwecken:
   if ~isempty(subroutine), feval(subroutine,x,n); end
   
   % Abbruchbedingung auswerten:
   fquer = sum(fx) / (n+1);
   abbruchbed = sqrt( sum( (fx-fquer).^2 )/n );
   
   iteration = iteration + 1;
%Ende der Iterationsschleife
end;

lsg = x(:,b);
