include <mach.h>
include <math/nlfit.h>
$if (datatype == r)
include "nlfitdefr.h"
$else
include "nlfitdefd.h"
$endif

define	COV	Mem$t[P2P($1)]	# element of COV


# NLERRORS -- Procedure to calculate the reduced chi-squared of the fit
# and the errors of the coefficients. First the variance
# and the reduced chi-squared of the fit are estimated. If these two
# quantities are identical the variance is used to scale the errors
# in the coefficients. The errors in the coefficients are proportional
# to the inverse diagonal elements of MATRIX.

procedure nlerrors$t (nl, z, zfit, w, npts, variance, chisqr, errors)

pointer	nl		# curve descriptor
PIXEL	z[ARB]		# data points
PIXEL	zfit[ARB]	# fitted data points
PIXEL	w[ARB]		# array of weights
int	npts		# number of points
PIXEL	variance	# variance of the fit
PIXEL	chisqr		# reduced chi-squared of fit (output)
PIXEL	errors[ARB]	# errors in coefficients (output)

int	i, n, nfree
pointer	sp, covptr
PIXEL   factor

begin
	# Allocate space for covariance vector.
	call smark (sp)
	call salloc (covptr, NL_NPARAMS(nl), TY_PIXEL)

	# Estimate the variance and reduce chi-squared of the fit.
	n = 0
	variance = PIXEL (0.0)
	chisqr = PIXEL (0.0)
	do i = 1, npts {
	    if (w[i] <= PIXEL (0.0))
	        next
	    factor = (z[i] - zfit[i]) ** 2
	    variance = variance + factor
	    chisqr = chisqr + factor * w[i]
	    n = n + 1
	}

	# Calculate the reduced chi-squared.
	nfree = n - NL_NFPARAMS(nl)
	if (nfree  > 0) {
	    variance = variance / nfree
	    chisqr = chisqr / nfree
	} else {
	    variance = PIXEL (0.0)
	    chisqr = PIXEL (0.0)
	}

	# If the variance equals the reduced chi_squared as in the case
	# of uniform weights then scale the errors in the coefficients
	# by the variance and not the reduced chi-squared

	if (abs (chisqr - variance) <= EPSILON$T) {
	    if (nfree > 0)
		factor = chisqr
	    else
		factor = PIXEL (0.0)
	} else
	    factor = 1.

	# Calculate the  errors in the coefficients.
	call aclr$t (errors, NL_NPARAMS(nl))
	call nlinvert$t (ALPHA(NL_ALPHA(nl)), CHOFAC(NL_CHOFAC(nl)),
	    COV(covptr), errors, PLIST(NL_PLIST(nl)), NL_NFPARAMS(nl), factor)

	call sfree (sp)
end


# NLINVERT -- Procedure to invert matrix and compute errors

procedure nlinvert$t (alpha, chofac, cov, errors, list, nfit, variance)

PIXEL	alpha[nfit,ARB]		# data matrix
PIXEL	chofac[nfit, ARB]	# cholesky factorization
PIXEL	cov[ARB]		# covariance vector
PIXEL	errors[ARB]		# computed errors
int	list[ARB]		# list of active parameters
int	nfit			# number of fitted parameters
PIXEL	variance		# variance of the fit

int	i, ier

begin
	# Factorize the data matrix to determine the errors.
	call nl_chfac$t (alpha, nfit, nfit, chofac, ier)

	# Estimate the errors.
	do i = 1, nfit {
	    call aclr$t (cov, nfit)
	    cov[i] = PIXEL (1.0)
	    call nl_chslv$t (chofac, nfit, nfit, cov, cov)
	    if (cov[i] >= PIXEL (0.0))
	        errors[list[i]] = sqrt (cov[i] * variance)
	    else
		errors[list[i]] = PIXEL (0.0)
	}
end