function [X,BestF,Iters] = hookejeeves(N, X, StepSize, MinStepSize, Eps_Fx, MaxIter, myFx)
% Function HOOKEJEEVS performs multivariate optimization using the
% Hooke-Jeeves search method.
%
% Input
%
% N - number of variables
% X - array of initial guesses 
% StepSize - array of search step sizes
% MinStepSize - array of minimum step sizes
% Eps_Fx - tolerance for difference in successive function values
% MaxIter - maximum number of iterations
% myFx - name of the optimized function
%
% Output
%
% X - array of optimized variables
% BestF - function value at optimum
% Iters - number of iterations
%

Xnew = X;
BestF = feval(myFx, Xnew, N);
LastBestF = 100 * BestF + 100;

bGoOn = true;
Iters = 0;

while bGoOn

  Iters = Iters + 1;
  if Iters > MaxIter
    break;
  end
  
  X = Xnew;
  
  for i=1:N
    bMoved(i) = 0;
    bGoOn2 = true;
    while bGoOn2
      xx = Xnew(i);
      Xnew(i) = xx + StepSize(i);
      F = feval(myFx, Xnew, N);
      if F < BestF
        BestF = F;
        bMoved(i) = 1;
      else
        Xnew(i) = xx - StepSize(i);
        F = feval(myFx, Xnew, N);
        if F < BestF
          BestF = F;
          bMoved(i) = 1;
        else
          Xnew(i) = xx;
          bGoOn2 = false;
        end  
      end  
    end
  end
  
  bMadeAnyMove = sum(bMoved);
  
  if bMadeAnyMove > 0 
    DeltaX = Xnew - X;
    lambda = 0.5;
    lambda = linsearch(X, N, lambda, DeltaX, myFx);
    Xnew = X + lambda * DeltaX;
  end
  
  BestF = feval(myFx, Xnew, N);
  
  % reduce the step size for the dimensions that had no moves
  for i=1:N
    if bMoved(i) == 0
      StepSize(i) = StepSize(i) / 2;
    end
  end
  
  if abs(BestF - LastBestF) < Eps_Fx
    break
  end
  
  LastBest = BestF;
  bStop = true;
  for i=1:N
    if StepSize(i) >= MinStepSize(i)
      bStop = false;
    end  
  end
  
  bGoOn = ~bStop;
  
end

function y = myFxEx(N, X, DeltaX, lambda, myFx)

  X = X + lambda * DeltaX;
  y = feval(myFx, X, N);

% end

function lambda = linsearch(X, N, lambda, D, myFx)

  MaxIt = 100;
  Toler = 0.000001;

  iter = 0;
  bGoOn = true;
  while bGoOn
    iter = iter + 1;
    if iter > MaxIt
      lambda = 0;
      break
    end  
   
    h = 0.01 * (1 + abs(lambda));
    f0 = myFxEx(N, X, D, lambda, myFx);
    fp = myFxEx(N, X, D, lambda+h, myFx);
    fm = myFxEx(N, X, D, lambda-h, myFx);
    deriv1 = (fp - fm) / 2 / h;
    deriv2 = (fp - 2 * f0 + fm) / h ^ 2;
    diff = deriv1 / deriv2;
    lambda = lambda - diff;
    if abs(diff) < Toler
      bGoOn = false;
    end
  end

% end

%{

Sub DirectSearch()
  Const MAX_VARS = 3
  
  Dim N As Integer, I As Integer
  Dim Row As Integer, MaxIter As Integer
  Dim X(MAX_VARS) As Double, Xnew(MAX_VARS) As Double
  Dim StepSize(MAX_VARS) As Double
  Dim MinStepSize(MAX_VARS) As Double
  Dim Deltax(MAX_VARS) As Double
  Dim F As Double, XX As Double, Lambda As Double
  Dim BestF As Double, LastBestF As Double, EPSF As Double
  Dim bStop As Boolean, bMadeAnyMove As Boolean
  Dim bMoved(MAX_VARS) As Boolean
  
  N = MAX_VARS
  For I = 1 To N
    Xnew(I) = Cells(I + 1, 2)
    StepSize(I) = Cells(I + 1, 3)
    MinStepSize(I) = Cells(I + 1, 4)
  Next I
  EPSF = Cells(2, 1)
  
  Range("E1:Z32000").Clear
  Range("E1").Value = "F(X())"
  Range("E1:Z1").Font.Bold = True
  For I = 1 To N
    Cells(1, 5 + I) = "X(" & I & ")"
    Cells(1, 5 + I + N) = "Step(" & I & ")"
  Next I
  ' calculate and display function value at initial point
  BestF = MyFx(Xnew, N)
  LastBestF = 100 * BestF + 100
  
  Row = 1
  Do
    Row = Row + 1
    
    For I = 1 To N
      X(I) = Xnew(I)
    Next I
    
    For I = 1 To N
      bMoved(I) = False
      Do
        XX = Xnew(I)
        Xnew(I) = XX + StepSize(I)
        F = MyFx(Xnew, N)
        If F < BestF Then
          BestF = F
          bMoved(I) = True
        Else
          Xnew(I) = XX - StepSize(I)
          F = MyFx(Xnew, N)
          If F < BestF Then
            BestF = F
            bMoved(I) = True
          Else
            Xnew(I) = XX
            Exit Do
          End If
        End If
      Loop
    Next I
    
    ' moved in any direction?
    bMadeAnyMove = True
    For I = 1 To N
      If Not bMoved(I) Then
        bMadeAnyMove = False
        Exit For
      End If
    Next I
    
    If bMadeAnyMove Then
      For I = 1 To N
        Deltax(I) = Xnew(I) - X(I)
      Next I
'''      Lambda = CubicIntMin(N, X, Deltax, 0, 1, 0.0001, 0.00001, 100)
'''      For I = 1 To N
'''        Xnew(I) = X(I) + Lambda * Deltax(I)
'''      Next I
      If LinSearch_DirectSearch(X, N, Lambda, Deltax, 0.1, 0.0001) Then
        For I = 1 To N
          Xnew(I) = X(I) + Lambda * Deltax(I)
        Next I
      End If
    End If
    BestF = MyFx(Xnew, N)
    
    ' reduce the step size for the dimensions that had no moves
    For I = 1 To N
      If Not bMoved(I) Then StepSize(I) = StepSize(I) / 2
    Next I
        
            
    ' show results
    Cells(Row, 5) = BestF
    For I = 1 To N
      Cells(Row, 5 + I) = X(I)
      Cells(Row, 5 + I + N) = StepSize(I)
    Next I
        
    ' test function value convergence
    'If Abs(BestF - LastBestF) < EPSF Then Exit Do
        
    LastBestF = BestF
    
    bStop = True
    For I = 1 To N
      If StepSize(I) >= MinStepSize(I) Then
        bStop = False
        Exit For
      End If
    Next I
    
  Loop Until bStop
  
End Sub


%}