function [Hssn,x0,f0,g0,dx,corr] = numhess(func,x0,points,delta,ix,diagon)
% [Hssn,x0,f0,dx,corr] = numhess(func,x0,{points,delta,ix,diagon}) 
%
% Numerically compute the Hessian of a function of several parameters.
% Use the Hessian to compute error bounds for the parameters.  
% Requires 1+2n^2 function evaluations, where n is the dimension of [ix].
% If the Hessian is not positive definite, re-run numhess() using the returned 
% value of x0.
%       -----  Input Values -------
% func   = function of parameters returning error criteria, J = func(x) 
% x0     = the optimal value of the parameters as found through minimization 
% points = the number of data points used in the minimization
% diagon = 1: compute only the diagonal elements of the Hessian, 0: all elements
% ix     = the elements of x0 for which statistics are to be determined
% delta  = fractional change in parameters for Hessian computation [0.02]
%       -----  Return Values -------
% Hssn   = Hessian matrix
% x0     = new optimal parameters, x, possibly better than the input x0
% f0     = new optimal value for f(x)
% dx     = asymptotic standard error in the parameters
% corr   = correlation matrix of the parameters
%
% H.P. Gavin, Civil & Environ. Eng'g, Duke Univ.,  2005-1-22

% default values 


if nargin < 3, points = length(x0)+1; end
if nargin < 4, delta = 0.02; end
if nargin < 5, ix= [1:length(x0)]; end
if nargin < 6, diagon = 0; end

 n = length(ix);
 x0 = x0(:);				% make column vector

 fprintf('initial value:\n');
 [f0, g0] = feval(func,x0)	

 m = length(g0);
 g = zeros(m,4);

 Hssn = zeros(n);			% initialize the Hessian

 for i=1:n
   for j=i:n
     if i==j				% compute diagonal elements only
        x(:,1:2) = [ x0 x0 ];   
	ii = ix(i);	jj = ix(j);
        delta_x = delta*(1+abs(x0(ii)));
        
        delta_x = delta;

        x(ii,1) = x0(ii) + delta_x;
        x(ii,2) = x0(ii) - delta_x;
        [f(1), g(:,1)] = feval(func,x(:,1));
        [f(2), g(:,2)] = feval(func,x(:,2));
        Hssn(i,i) = (f(1) - 2*f0 + f(2)) / delta_x^2;

	[f,idx] = sort(f(1:2));
	x = x(:,idx);
	g = g(:,idx);       
	if (f(1) < 0.98*f0) & (max(g(:,1)) < 0)
		f0 = f(1);
		x0 = x(:,1);
 		fprintf('**** optimum updated ... f0 = %f  max_g0 = %f \n', ...
                                          f0, max(g0) );
	end

     elseif (diagon==0)			% compute off-diagonal elements only
	x(:,1:4)= [ x0 x0 x0 x0 ]; 

	ii = ix(i);	jj = ix(j);

        delta_xi = delta*(1+abs(x0(ii)));
        delta_xj = delta*(1+abs(x0(jj)));

        delta_xi = delta;
        delta_xj = delta;

	x(ii,1) = x0(ii) + delta_xi;
	x(ii,2) = x0(ii) + delta_xi;
	x(ii,3) = x0(ii) - delta_xi;
	x(ii,4) = x0(ii) - delta_xi;

	x(jj,1) = x0(jj) + delta_xj;
	x(jj,2) = x0(jj) - delta_xj;
	x(jj,3) = x0(jj) + delta_xj;
	x(jj,4) = x0(jj) - delta_xj;

        [f(1), g(:,1)] = feval(func,x(:,1));
        [f(2), g(:,2)] = feval(func,x(:,2));
        [f(3), g(:,3)] = feval(func,x(:,3));
        [f(4), g(:,4)] = feval(func,x(:,4));

	Hssn(i,j) = ( f(1) - f(2) - f(3) + f(4) ) / (4*delta_xi*delta_xj);
	Hssn(j,i) = Hssn(i,j);

	[f,idx] = sort(f(1:4));
	x = x(:,idx);
	g = g(:,idx);       
	if (f(1) < 0.98*f0) & (max(g(:,1)) < 0)
		f0 = f(1);
		x0 = x(:,1);
	        g0 = g(:,1);
 		fprintf('**** optimum updated ... f0 = %f  max_g0 = %f \n', f0, max(g0) );
	end

     end
   end
 end

 if any( eig(Hssn) < 0 ) 		% Hssn should be positive definite
	disp('  warning: This Hessian matrix is not positive definite.');
	disp('           If this is not a constrained solution,');
	disp('           re-start numhess about the new minimum point:');
 end

 measured_noise = sqrt(abs(f0)/(points-n+1));		% sqrt(SSR/dof)
 covar = 2*inv(Hssn) * abs(f0)/(points-n+1);   % f = sum((y_data - y_fit).^2);

 dx = sqrt(diag(covar)); 

 corr = covar ./ [dx * dx'];

 format short e
 disp('param no.     param            asymptotic standard error ')
 disp('=========     ==========    =============================== ')
 for i=1:length(ix)
     fprintf('   %2d        %12.5e   %12.5e   (+/- %6.1f percent)\n', ...
                              ix(i), x0(ix(i)), dx(i), 100*dx(i)/x0(ix(i)) );
 end

% endfunction # ------------------------------------------------------ NUMHESS