% J. C. Spall, March 2000
% This code computes the FDSA estimate for a given no. of loss measurements.  
% Code allows multiple
% replications ("cases") for averaging solutions.
% Separate code is given in Chapter 7, which is designed for Monte Carlo comparisons of 
% FDSA and SPSA.
% This code includes a d vector for use in storing data for plotting or analysis purposes. 
% This vector has nothing per se to do with the algorithm or optimization process (it was
% used in creating Figure 6.2 in ISSO).
%
global sigma p noise
p=2;
nFD=100;			% No. of loss measurements used in FDSA (per realization)
sigma=1;
loss='loss_waste_noise';
loss_nonoise='loss_waste_noisefree';
cases=50; 
alpha=1;
gamma =.16666701;
% Gain coefficients 
aFD=.4;
cFD=1;
AFD=0;
e=eye(p);
% INITIALIZATION
theta_0=3*ones(p,1); 
theta_lo=0*ones(p,1);   	%lower bounds on theta  
theta_hi=5*ones(p,1);    	%upper bounds on theta 
rand('seed',31415927)
randn('seed',511111113)
meanlossFD=0;
errthetaFD=0;
losssqFD=0;
ghat=zeros(p,1);		%dummy statements for setting dimension
yplus=zeros(1,1);
yminus=zeros(1,1);
t=2.0096;				%t-quantile for confidence interval calculation
truetheta=[1,0.666667]';
for i=1:cases
  replication=i 	% for monitoring output in large-scale runs 
  theta=theta_0;
  d(1,i,2)=feval(loss_nonoise,theta); 	%true loss value for evaluation and plotting purposes 
  													%(not part of optimization process)
  for k=1:nFD/(2*p)
    ak = aFD/(k+AFD)^alpha;
    ck = cFD/k^gamma;
    for j=1:p
      thetaplus = theta+ck*e(:,j);
      thetaminus = theta-ck*e(:,j);
      yplus=feval(loss,thetaplus);
      yminus=feval(loss,thetaminus);
      ghat(j)=(yplus - yminus)/(2*ck);
    end
    theta=theta-ak*ghat;
    % Two lines below invoke constraints
    theta=min(theta,theta_hi);
    theta=max(theta,theta_lo);
    d(k+1,i,2)=feval(loss_nonoise,theta); %records true loss value for ith case, kth iteration  
  end
  eval=feval(loss_nonoise,theta);
  meanlossFD=(i-1)*meanlossFD/i+eval/i;
  losssqFD=(i-1)*losssqFD/i+eval^2/i;
  errthetaFD=(i-1)*errthetaFD/i+((theta-truetheta)'*(theta-truetheta))^.5/i;  
end
meanlossFD
errthetaFD
if cases > 1
   sd_FD=(cases^(-.5))*((cases/(cases-1))^.5)*(losssqFD-meanlossFD^2)^.5
   % Confidence intervals
   [meanlossFD-t*sd_FD,meanlossFD+t*sd_FD]
else   
end
