% Haerer, Bernhardt, Schulz and Corripio (2012)
% "PRACTISE - Photo Rectification And ClassificaTIon SoftwarE"
%
%   written by
%   Karsten Schulz, Stefan Härer (LMU Munich)
%   08/2012
%   contact: s.haerer@iggf.geo.uni-muenchen.de
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Xopt, Fopt] = DDS_PRACTISE(FUN, X, LB, UB, R, MaxEval, fout_folder, dem, header)
%   Name:       DDS_PRACTISE
%   Purpose:    Dynamically Dimensioned Search (DDS) Algorithms searching 
%               for a global minimum of any objective function (Tolson & 
%               Shoemaker, 2007) adapted to PRACTISE
%
%   Output:     Fopt (Value of the objective function in the optimum)
%               Xopt (Parameters for optimal objective function)
%   Input:      FUN (Handle to the objective Function)
%               X (Initial guess for the optimal parameters)
%               LB (Lower boundary of the parameter space)
%               UB (Upper boundary of the parameter space)
%               R (Neighborhood pertubation size)
%               MaxEval (Maximum number of function evaluations)
%               fout_folder (Folder directory for DDS optimization path figure)
%               dem (DEM raster from ASCII-file [cf. ESRI ArcGIS])
%               header (values od DEM ASCII-file header [cf. ESRI ArcGIS]:
%                 ncols, nrows, xllcorner, yllcorner, cellsize, NODATA_value)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1: DDS parameters
%   specified in parameter file
%
% Step 2: initial evaluation of the function
disp('The number of DDS runs evaluated is')
count   = 1;
[F, gcpP] = feval(FUN, X, dem, header); 
%   plot initial projected GCP 
plot(gcpP(1,:)+0.5, gcpP(2,:)+0.5, 'o', 'MarkerEdgeColor', [1 0 0], 'MarkerSize',8)
drawnow
%
Fopt    = F;
Xopt    = X;
quit    = 0;
%
while(quit==0)
%    
    % Step 3: Randomly select J of the N variables
    Pi = 1-log(count)/log(MaxEval);
    Pn = rand(size(X));
    J=find(Pn<Pi);
    nj=length(j);
    if (nj==0)
        J=floor(rand(1)*N+1);
    end
    %
    % Step 4: calculate updated xnew
    sig    = R*(UB-LB).*randn(size(X));
    Xnew   = Xopt;
    Xnew(J)= Xnew(J)+sig(J);
    %
    i      = find(Xnew<LB);
    Xnew(i)= LB(i)+(LB(i)-Xnew(i));
    i1     = find(Xnew(i)>UB(i));
    Xnew(i1)=LB(i1);
    %
    z       = find(Xnew>UB);
    Xnew(z) = UB(z) - (Xnew(z)-UB(z));
    z1      = find(Xnew(z)<LB(z));
    Xnew(z1)= UB(z1);
    %
    % Step 5: 
    count   = count+1;
    [F,gcpP]       = feval(FUN, Xnew, dem, header);
    if (F<Fopt)
            Fopt    = F
            Xopt    = Xnew;
%           % plot improved projected GCPs
%           plot(gcpP(1,:)+0.5, gcpP(2,:)+0.5, '.', 'MarkerEdgeColor', [1 0 0], 'MarkerSize',4)
%           drawnow
    end
    %
    % Step 6: Stop Criteria
    if(mod(count,500)==0)
        disp(count)
    end    
    if(count==MaxEval)
        quit = 1;
    end    
%
% end of while loop
%
end




