function ddEUU = lyngby_nn_qddeus(X, V, W, H, H12)

% lyngby_nn_qddeus     - 2nd der., quad., all weights, gauss approx.
%
%       function ddEUU = lyngby_nn_qddeus(X, V, W, H, H12)
%
%	Input:  X       Neural network input
%               V       Input weights
%               W       Output weights
%               H12     (1-tanh(XV).^2) 
%
%       Output: ddEUU   Derivative as a Hessian matrix
%
%       Gauss-Newton (symmetric) approximation to the second order
%       derivative of the qaudratic error function with respect to
%       both the input and the output weights in the neural network.

% cvs : $Id: lyngby_nn_qddeus.m,v 1.1 1998/02/13 16:59:29 fnielsen Exp $
%       $Revision: 1.1 $


    [Np, Nii] = size(X);
    [Nii, Nh] = size(V);
    [Nhh, No] = size(W);
    Nvw1 = Nii * Nh;
    Nvw2 = Nhh * No;
    Nvw = Nvw1 * Nvw2;

    if Nhh > Nh
      H = [H ones(Np,1)];
    end
    
    if No == 1
      w = W(1:Nh, :)';
    else
      w = sum(W(1:Nh, :)');
    end

    A = kron(ones(1,Nh), X) .* kron(H12 .* (ones(Np,1) * w), ones(1, Nii));
    ddEVV = A'*A / (Np*No);
    ddEWW = lyngby_nn_qddews(H, W);
    
    indexW1 = fix((0:(Nvw1-1)) / Nii)+1;
    indexW2 = fix((0:(Nvw2-1)) / Nhh)+1;
    W1 = W(indexW1,:)'; 
    W1 = W1(indexW2,:); 
    w = W1(:)';
    indexX   = rem(fix((0:(Nvw-1)) / (Nhh*No)), Nii) + 1;
    indexH12 = fix((0:(Nvw-1)) / (Nii*Nhh*No)) + 1;
    indexH   = rem((0:(Nvw-1)), Nhh)+1;

    A = X(:, indexX) .* H12(:, indexH12) .* (ones(Np, 1) * w) .* ...
	H(:, indexH);

    ddEVW = reshape(sum(A) / (No*Np), Nvw2, Nvw1);
    
    ddEUU = [ ddEVV ddEVW' ; ddEVW ddEWW ];
    
    if 1==0
      ddEVVtest = zeros(Nii*Nh);
      for p=1:Np
	for o=1:No
	  for h1=1:Nh
	    for h2 = 1:Nh
	      for i1 = 1:Nii
		for i2 = 1:Nii
		  n1 = i1+Nii*(h1-1);
		  n2 = i2+Nii*(h2-1);
		  ddEVVtest(n1,n2) = ddEVVtest(n1,n2) + ...
		      1/(Np*No) * W(h1,o)*W(h2,o) * H12(p,h1)*H12(p,h2) ...
		      * X(p,i1)*X(p,i2);
		end
	      end
	    end
	  end
	end
      end
    
      ddEVWtest = zeros(Nii*Nh, Nhh*No);
      for p=1:Np
	for o1=1:No
	  for h1=1:Nhh
	    for h2 = 1:Nh
	      for i2 = 1:Nii
		n1 = i2+Nii*(h2-1);
		n2 = h1+Nhh*(o1-1);
		ddEVWtest(n1,n2) = ddEVWtest(n1,n2) + ...
		    1/(Np*No) * H(p,h1) * W(h2,o1) * H12(p,h2) ...
		    * X(p,i2);
	      end
	    end
	  end
	end
      end

      ddEUUtest = [ ddEVVtest ddEVWtest ; ddEVWtest' ddEWW ];

    end






