function [g,ak,delta,ck]= SPSA(func,x,nrep,nci,a,c,maxiter,niter)
%SPSA average gradient for nrep individual gradients
%nrep is the number of individual gradients used to get avg grad
% x is the current iterate at which grad is to be calculated
%nci is the dimenion of the problem
% a is the parameter used to calculate step size
% c is parameter used to calculate perturbation
tic
xlow=zeros(nci,1);
xup=ones(nci,1);
%SPSA Parameters (typical values used)
alpha2=0.602;
gamma=0.101;
A=0.1*maxiter;
% a=0.16;
% c=0.1;

sum=0;

ak=a/(niter+1+A)^alpha2;
ck=c/(niter+1)^gamma;

for i=1:nrep
    delta(:,i)=2*round(rand(nci,1))-1;
    xplus(:,i)=x+ck*delta(:,i);
    xminus(:,i)=x-ck*delta(:,i);

end

parfor i=1:nrep
    yplus(i)=func(xplus(:,i)',i);
    yminus(i)=func(xminus(:,i)',i);
    
end
for i=1:nrep
    g(:,i)=(yplus(i)-yminus(i))./ (2*ck*delta(:,i));
    sum=sum+g(:,i);
end
g=sum/nrep;

toc
end

