% J. C. Spall, March 1999
% This code provides values of a, A, and c in the gains a_k=a/(k+1+A)^alpha and 
% c_k = c/(k+1)^.101.  Code assumes Bernoulli +/- 1 distribution for the delta elements. 
%
% Note: if get some error "Subscript indices must be integer values"
% then use "clear" function and redo.
%
% Specify the dimension, loss function, and i.c. for theta; also alpha for gain 
% selection (alpha = 0.602 is often used as a companion to gamma = 0.101). 
%
global sigma p z	% sigma is noise level used in loss calls; z is global 
						% random variable for use in generating noise in loss
						% function calls(so that seed can be set here).
sigma=.1;
p=10;
loss='loss_test';
theta=6*[1,1,1,1,1,1,1,1,1,1]';
alpha=0.602;
%
% User input on measurement noise level, expected no. of iterations in the SPSA run,
% desired magnitude of change in the theta elements, the number of SPSA gradient approximations
% that will be averaged, and the no. of loss evaluations to be used in the gain calculations
% here (note this no. should be divisible by twice the no. of averaged gradients).
%
step= input('What is the initial desired magnitude of change in the theta elements? ');
gavg= input('How many averaged SP gradients will be used per iteration? ');
A = .10*input('What is the expected number of loss evaluations per run? ')...
   /(2*gavg);
c = input('What is the standard deviation of the measurement noise at i.c.? ');
c = max(c/gavg^0.5, .0001); %averaging set up for independent noise
NL = input('How many loss function evaluations do you want to use in this gain calculation? ');
%
% Do the NL/(2*gavg) SP gradient estimates
%
rand('seed',31415927)	% Seed for delta generation
randn('seed',111113); 	% Seed for noise in loss measurements
gbar=zeros(p,1);
for i=1:NL/(2*gavg)
   ghat=zeros(p,1);
   for j=1:gavg
      delta=2*round(rand(p,1))-1;
      thetaplus=theta+c*delta;
      thetaminus=theta-c*delta;
      yplus=feval(loss,thetaplus);
      yminus=feval(loss,thetaminus);
      ghat=(yplus-yminus)./(2*c*delta)+ghat;
   end
   gbar=gbar+abs(ghat/gavg);
end
gbar
meangbar=mean(gbar);
meangbar=meangbar/(NL/(2*gavg));
a=step*((A+1)^alpha)/meangbar;
%
% Print out of c, A, and a
%
c
A
a