%Version 1.0 Dynamic Function #1 (DF1)
% Ronald W. Morrison, June 1999
%
% This is MATLAB code that generates a "field of moving cones"
% for use as a dynamic test function for testing evolutionary algorithms.
% It is described in the paper "A Test Problem Generator for 
% Non-Stationary Enviornments" by Ron Morrison and Ken De Jong
% in the Proceedings of the 1999 Congress on 
% Evolutionary Computation (CEC-99).  There is also a small user's
% guide available with this code at www.aic.nrl.navy.mil

clear;

% Set the number of peaks,N
% Rbase sets the basic "spikiness" of the peaks, big = spiky,
% Rrange sets the range of the spikiness among the peaks, 
% Hbase sets the basic height of the peaks,
% Hrange sets the range of the difference is height off of the base,
% Notation convention: X,Y are the plotting matrices, 
%                      x,y are the origin parameter matrices
%
%set these to whatever values you would like
N = 15;
Rbase = 8;
Rrange = 12;
Hbase = 1;
Hrange = 9;
%
%
%the max & mins are computed here 
Rmax = Rbase + Rrange;
Rminpct = Rbase/Rmax;
Hmax = Hbase + Hrange;
Hminpct = Hbase/Hmax;

%This program will make a MATLAB movie of the movement.
% Numframes sets the number of iterations and also the
% number of frames in the movie,

NumFrames = 30;

M = moviein(NumFrames);

% The dynamics can be in height, "spikiness", x location, Y location,
% or any combincation of these dynamics.

% Now select whether the movement is in Height (H), Spikiness(R), 
% X location (x), Y location (y), or any combincation.  
% Set the ones you want to "1". (for example, the settings shown
% here will provide dynamics in x and y only).

Hboolean = 0;
Rboolean = 0;
xboolean = 1;
yboolean = 1;

% Not all of the peaks need move. 
% Now set NumMove, the number of peaks that will move (NumMove)
% to a number less than or equal to N.

NumMove = 14;

% Now set A's for motion (different for each
%     characteristic); 
% A is the constant for the Logistics Function (see paper).
% There is also a scaling factor in case the
% step values that the Logistics Function generates are
% not in the range that you like.
% (0 < step scale factor <.5)

%examples: 

%for height:
Ah = 1.2;
Hstepscale = .3;
% for R:
Ar = 1.6;
Rstepscale = .5;
% for x:
Ax = 3.8;
xstepscale = .5;
% for y:
Ay = 3.8;
ystepscale = .5;


% Now select the cones to apply the motion to with a vector C of indices from 1 to N


% C indices are the specific cones that you want to move.  For 
% example if you were only moving 2 cones and you wanted to move
% cones 4 and 12, you would put C = [4 12].  In general this is
% of limited usefulness, just MAKE SURE that you have as many 
% indices as you have moving cones.
 
C = [2 3 4 5 6 7 8 9 10 11 12 13 14 15];

% Now, if you want, you can set some cones to start their movement by getting larger
% and some by getting smaller.  1 means start by getting larger.  0 means
% start by getting smaller.  AGAIN, YOU MUST HAVE ONE VALUE FOR EACH MOVING CONE>

% IUP sets the initial movement as up or down, with reversal at the max, for 
% each parameter H, R, x, or y.  

IUPH = [0 1 0 1 0 0 1 1 1 0 1 0 1 0]; 
IUPR = [0 1 0 1 0 0 1 1 1 0 1 0 1 0]; 
IUPx = [1 1 1 1 1 1 1 1 1 1 1 1 1 1]; 
IUPy = [1 0 1 0 1 0 1 0 1 0 1 0 1 0]; 

% The logistics function has some settling required for the initial step value.
% That is done next.

%Now work past the initial transients in the step values
  Hstep = .5;
  for i=2,100
     Hstep = Ah*Hstep*(1-Hstep);
  end
  Rstep = .5;
  for i=2,100
     Rstep = Ar*Rstep*(1-Rstep);
  end
  xstep = .5;
  for i=2,100
     xstep = Ax*xstep*(1-xstep);
  end
  ystep = .5;
  for i=2,100
     ystep = Ay*ystep*(1-ystep);
  end
  
%The landscape is between -1 and 1 on the XY plane
% establish the X-Y range.  The plotting matrix is XY.
% The picture is drawn with the meshgrid Matlab function.

[X,Y] = meshgrid(-1:.01:1);

%initialize radius, height, and x(0) and y(0) for all the cones
% You can create specific peaks here, just hard code them
% and adjust the loop accordingly.

% Hard code cone #1 to al floor of H(1)=0 and R(1)= 0 if you do NOT want
% negative fitnesses. You can also use this technique to select a 
% minimum other than 0. 

R(1) = 0;
H(1) = 0;
x(1) = .5;
y(1) = .5;

% randomly initialize the rest (unless you hard coded some others. In that
% case change the indes so as to not re-randomize cones that you hard coded.

for i = 2:N
  R(i) = Rbase + Rrange*rand;
  H(i) = Hbase + Hrange*rand;
  x(i) = 2*rand - 1;
  y(i) = 2*rand - 1;
end

%This plots the first cone.
 
Z =(H(1) - R(1)*(sqrt((X-x(1)).^2 + (Y-y(1)).^2)));

% This loop successively takes the MAX of all of the cones to build a surface.

for i = 2:N
  ZI =(H(i) - R(i)*(sqrt((X-x(i)).^2 + (Y-y(i)).^2)));
  Z = max(Z,ZI);  
end

% This draws the picture
meshc(Z);
axis([0 220 0 220 0 15]); %The Z axis is set to the maximum height!

%now apply the motion (height, R, x, or y)

for m = 1:NumFrames

%first move in H ... if you have hboolean set.
   if (Hboolean == 1)
   for L = 1:NumMove;
     Hpct = H(C(L))/Hmax;               % convert H value to a percent of max 
     Hstep = Ah*Hstep*(1-Hstep);
     Htemp = Hstep*Hstepscale;          %scale the step
     if (IUPH(L)==1)                     % if on the increase
       Hpct = Hpct + Htemp;            % increase by chosen A and step size 
       if (Hpct >1.0)                   % too large -- reverse
         IUPH(L) = 0;
         Hpct = Hpct - 2*Htemp;
       end
     else                                %decrease by chosen A and step size
       Hpct = Hpct - Htemp;
       if (Hpct < Hminpct)                    % too small, reverse
         IUPH(L) = 1;        
         Hpct = Hpct + 2*Htemp;
       end
     end
     H(C(L)) = Hpct * Hmax;
   end
   end   % moving H  
% Now move in R
   if (Rboolean == 1)
   for L = 1:NumMove;
     Rpct = R(C(L))/Rmax;               % convert H value to a percent of max 
     Rstep = Ar*Rstep*(1-Rstep);   
     Rtemp = Rstep * Rstepscale;    %scle the step
     if (IUPR(L)==1)
       Rpct = Rpct + Rtemp;            % increase by chosen A and step size 
       if (Rpct >1.0)                   % too large -- reverse
         IUPR(L) = 0;
         Rpct = Rpct - 2*Rtemp;
       end
     else                                %decrease by chosen A and step size
       Rpct = Rpct - Rtemp;
       if (Rpct < Rminpct)                    % too small, reverse
         IUPR(L) = 1;        
         Rpct = Rpct + 2*Rtemp;
       end
     end
     R(C(L)) = Rpct * Rmax;
   end
   end % moving R
% Now move in x
   if (xboolean == 1)
   for L = 1:NumMove;
     xpct = (x(C(L))+1)/2;               % convert X value to a percent of max 
     xstep = Ax*xstep*(1-xstep);
     xtemp = xstep * xstepscale;   %scale the step     
     if (IUPx(L)==1)
       xpct = xpct + xtemp;            % increase by chosen A and step size 
       if (xpct >1.0)                   % too large -- reverse
         IUPx(L) = 0;
         xpct = xpct - 2*xtemp;
       end
     else                                %decrease by chosen A and step size
       xpct = xpct - xtemp;
       if (xpct < 0)                    % too small, reverse
         IUPx(L) = 1;        
         xpct = xpct + 2*xtemp;
       end
     end
     x(C(L)) = (xpct * 2) - 1;
   end
   end   % moving x
% Now move in y
   if (yboolean == 1) 
   for L = 1:NumMove;
     ypct = (y(C(L))+1)/2;               % convert X value to a percent of max 
     ystep = Ay*ystep*(1-ystep);
     ytemp = ystep * ystepscale;         %scale the step     
     if (IUPy(L)==1)
       ypct = ypct + ytemp;            % increase by chosen A and step size 
       if (ypct >1.0)                   % too large -- reverse
         IUPy(L) = 0;
         ypct = ypct - 2*ytemp;
       end
     else                                %decrease by chosen A and step size
       ypct = ypct - ytemp;
       if (ypct < 0)                    % too small, reverse
         IUPy(L) = 1;        
         ypct = ypct + 2*ytemp;
       end
     end
     y(C(L)) = (ypct * 2) - 1;
   end
   end   % moving y

%now plot the revised cones   
       Z = (H(1) - R(1)*(sqrt((X-x(1)).^2 + (Y-y(1)).^2)));
   for i = 2:N
       ZI =(H(i) - R(i)*(sqrt((X-x(i)).^2 + (Y-y(i)).^2)));
       Z = max(Z,ZI);  
   end
   meshc(Z);
   axis([0 220 0 220 0 15]); %The Z axis is the maximum height
   M(:,m) = getframe;        %This makes a Matlab movie if you want one
  %pause                     %This pauses Matlab if you want to study the shape
end