%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kevin Passino
%   Version: 4/2/01
% 
% This program simulates the minimization of a simple function with
% the SPSA method.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear                % Initialize memory

p=2;                 % Dimension of the search space 

thetamin=[0; 0];  % Set edges of region want to search in
thetamax=[30;30];


Nspsa=100;            % Maximum number of iterations to produce

% Next set the parameters of the algorithm:

lambda=1;
lambda0=1;

alpha1=0.602;
alpha2=0.101;

c=0.01;
		  
% Next, pick the initial value of the parameter vector

theta(:,1)=[15; 15];  % Falls into one local minimum
%theta(:,1)=[17; 10];  % Moves to the global minimum
%theta(:,1)=[15; 20];  % Start at global maximum - goes to a local minimum
%theta(:,1)=[1; 1];  % Start at a relatively flat region - moves slow


% Allocate memory 

theta(:,2:Nspsa)=0*ones(p,Nspsa-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start the SPSA loop

for j=1:Nspsa

	% Use projection in case update (or initial values) out of range
	
	theta(:,j)=min(theta(:,j),thetamax);
	theta(:,j)=max(theta(:,j),thetamin);

	% Set parameters and perturb parameters
	
	lambdaj=lambda/(lambda0+j)^alpha1;  % Note since this corresponds to zero
	cj=c/j^alpha2;
	Delta=2*round(rand(p,1))-1; % According to a Bernoulli +- 1 dist.
	thetaplus=theta(:,j)+cj*Delta;
	thetaminus=theta(:,j)-cj*Delta;
	
	% Use projection in case perturbed out of range.
	
	thetaplus=min(thetaplus,thetamax);
	thetaplus=max(thetaplus,thetamin);
	thetaminus=min(thetaminus,thetamax);
	thetaminus=max(thetaminus,thetamin);

	% Next we must compute the cost function for thetaplus and thetaminus
	
	Jplus=optexampfunction([thetaplus(1,1);thetaplus(2,1)]);
	Jminus=optexampfunction([thetaminus(1,1);thetaminus(2,1)]);
		
	% Next, compute the approximation to the gradient
	
	g=(Jplus-Jminus)./(2*cj*Delta);
	
	% Then, update the parameters
	
	theta(:,j+1)=theta(:,j)-lambdaj*g;
		
end % End main loop...




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the function we are seeking the minimum of:

x=0:31/100:30;   % For our function the range of values we are considering
y=x;

% Compute the function that we are trying to find the minimum of.

for jj=1:length(x)
	for ii=1:length(y)
		z(ii,jj)=optexampfunction([x(jj);y(ii)]);
	end
end

% First, show the actual function to be maximized

figure(1)
clf
surf(x,y,z);
colormap(jet)
% Use next line for generating plots to put in black and white documents.
colormap(white);
xlabel('x=\theta_1');
ylabel('y=\theta_2');
zlabel('z=J');
title('Function to be minimized');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, provide some plots of the results of the simulation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t=0:Nspsa;  % For use in plotting

figure(2) 
clf
plot(t,theta(1,:),'k-',t,theta(2,:),'k--')
ylabel('\theta_1, \theta_2')
xlabel('Iteration, j')
title('SPSA parameter trajectories')


figure(3) 
clf
contour(x,y,z,25)
colormap(jet)
% Use next line for generating plots to put in black and white documents.
colormap(gray);
xlabel('x=\theta_1');
ylabel('y=\theta_2');
title('Function to be minimized (contour map)');

hold on

plot(theta(1,:),theta(2,:),'k-')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%