%function for computing gradient of theta 
% z = (nu,mu,u1,u2,u3,u4,psi,x)
function [gradF,gradG] = gradfun(z)
z = z(:);
nu = z(1:4);
mu = z(5:7);
U = [z(8:11) z(12:15) z(16:19) z(20:23)];
x = z(24:25);
psimax = z(26);
y = 0;
  for j = 1:4,
    y = y + nu(j)*gradgx(x,U(:,j));
  end %for
y = y + mu(1)*[x(2);x(1)] + mu(2)*[1;-2] + mu(3)*[-1;0.5];
%derivatrives with respect to nu
for i = 1:4,
  gradFnu(i) = psimax - g(x,U(:,i)) + y'*gradgx(x,U(:,i));
%derivatives with respect to ui
  u = U(:,i);
  v(1) = 400*u(2) + 2000;
  v(2) = 250*u(3) + 1000;
  v(3) = 100*u(1) + 500;
  v(4) = exp(0.7386 + 0.0998*u(4));
  Dv_Du = [0 400 0 0
           0  0  250 0
           100 0 0 0
           0 0 0 0.0998*exp(0.7386 + 0.0998*u(4))];
  D2Gxv(1,1) = 4/((x(1)*x(2))^2*v(4));
  D2Gxv(2,1) = 8/(x(1)^3*x(2)*v(4));
  D2Gxv(3,1) = 4*v(3)/(x(1)^3*(x(2)*v(4))^2);
  D2Gxv(4,1) = -4*v(1)/((x(1)*x(2)*v(4))^2)...
               -8*v(2)/(x(1)^3*x(2)*v(4)^2)...
               -4*v(3)^2/((x(1)*v(4))^3*x(2)^2);
  D2Gxv(1,2) = 8/(x(1)*x(2)^3*v(4));
  D2Gxv(2,2) = 4/((x(1)*x(2))^2*v(4));
  D2Gxv(3,2) = 4*v(3)/(x(2)^3*(x(1)*v(4))^2);
  D2Gxv(4,2) = -8*v(1)/(x(1)*x(2)^3*v(4)^2)...
               -4*v(2)/((x(1)*x(2))^2 *v(4)^2)...
               -4*v(3)^2/((x(2)*v(4))^3*x(1)^2);
% signs reversed to account that we are dealing with g_tilde = -g
  D2Gxu = -Dv_Du'*D2Gxv;
  gradFu(:,i) = nu(i)*(- gradgu(x,U(:,i)) + D2Gxu*y);
end %for
gradFnu = gradFnu(:);
gradFu = gradFu;
gradFmu(1) = -(x(1)*x(2) - 300) + y'*[x(2);x(1)];
gradFmu(2) = -(x(1) - 2*x(2)) + y'*[1;-2];
gradFmu(3) = -(0.5*x(2) - x(1)) + y'*[-1;0.5];
gradFmu = gradFmu(:);
gradF = [gradFnu;gradFmu;gradFu(:,1);gradFu(:,2);gradFu(:,3);gradFu(:,4);0;0;0];
l = ones(7,1);
for i=1:4,
  G = 2*U(:,i);
  gradG(:,i) = [zeros(7+3*(i-1),1);G;zeros(15-3*(i-1),1)];
end %for
gradG(:,5) = [l;zeros(19,1)];
gradG(:,6) = [-l;zeros(19,1)];