% 3.
maxHiddenNeuronCount = 10;

for hiddenNeuronCount = 1:maxHiddenNeuronCount
   net = train_net(trainData(:,1), trainData(:,2), hiddenNeuronCount);
   Z = calc_jacobian(net, trainData);
   rankZ(hiddenNeuronCount) = rank(Z);
   q(hiddenNeuronCount) = 2*net.layerweights{2,1}.size(2)+(net.layerweights{2,1}.size(2)+1)*1;

   [U W V] = svd(Z);
   VW = V./repmat([diag(W)]', size(V,1), 1);
   ZTZ1 = inv(VW*VW)*V';
   diagH(:, hiddenNeuronCount) = diag(Z * ZTZ1 * Z');
   hVariance(hiddenNeuronCount) = var(diagH(:, hiddenNeuronCount));

   if(q(hiddenNeuronCount) ~= rankZ(hiddenNeuronCount))
       hVariance(hiddenNeuronCount) = 0;
   end

   if(hVariance(hiddenNeuronCount) > 1)
       hVariance(hiddenNeuronCount) = 1;
   end
end