function ddEV = lyngby_nn_eddevds(X, Y, V, W, H12)

% lyngby_nn_eddevds    - 2nd order, entropic, input, diag. sym.
%
%       function ddEV = lyngby_nn_eddevds(X, Y, V, W, H12)
%
%       Input:  X      Neural network input
%               Y      Neural network output
%               W      Output weights
%               H12    One minus Hidden layer in 2nd power
%
%       Output: ddEV   Derivative, a vector
%
%       'Symmetric' diagonal approximation to the second order
%       derivative of the entropic error function with respect 
%       to the input weights.
%
%       H12 can be computed as (1-H.^2)
%
%       See also: lyngby_nn_eddewd

% cvs : $Id: lyngby_nn_eddevds.m,v 1.1 1997/08/22 16:56:17 fnielsen Exp $
%       $Revision: 1.1 $

    [Nii, Nh] = size(V);
    [Np, No] = size(Y);

    H122  = H12(:,1:Nh) .* H12(:,1:Nh);  % Np x Nh
    X2t   = (X .* X)';                   % Nii x Np
    Y12   = 1 - Y.*Y;                    % Np x No
    W2t   = (W(1:Nh,:) .* W(1:Nh,:))';   % No x Nh
    
    ddEV  = 1/(No*Np) * X2t * ((Y12 * W2t) .* H122);
    
    % Y12W2 = Y12 * W2t;               % Np x Nh
    % WYH   = Y12W2 .* H122(1:Nh,:);   % Np x Nh 
    % ddEV  = 1/(No*Np) * X2t * WYH;   % Nii x Nh
    
    % Numerical errors? 
    indices = find(V==0);
    if indices
      ddEV(indices) = zeros(1,length(indices));
    end

    % This is equivalent to the last computation:
    % ddEV = zeros(Nh, Nii);
    % for oo = 1:No
    %   for hh = 1:Nh
    %     for pp = 1:Np
    %       ddEV(hh,:) = W2(oo,hh) * Y12(oo,pp) * H122(hh,pp) * X2t(pp,:);
    %	  end
    %   end
    % end















