Initial commit

1 files changed, 121 insertions, 0 deletions

diff --git a/pkg/xtools/inlfit/indeviantd.x b/pkg/xtools/inlfit/indeviantd.x
new file mode 100644
index 00000000..ec32e637
--- /dev/null
+++ b/pkg/xtools/inlfit/indeviantd.x
@@ -0,0 +1,121 @@
+include	<mach.h>
+
+
+# IN_DEVIANT -- Find deviant points with large residuals from the fit
+# and reject them from the fit.  The sigma of the fit residuals is calculated.
+# The rejection thresholds are set at (+/-)reject*sigma.  Points outside the
+# rejection threshold are recorded in the reject array.
+
+procedure in_deviantd (nl, x, y, w, rejpts, npts, nvars, low_reject,
+	high_reject, grow, nreject, newreject)
+
+pointer	nl				# NLFIT descriptor
+double	x[ARB]				# Input ordinates (npts * nvars)
+double	y[npts]				# Input data values
+double	w[npts]				# Weights
+int	rejpts[npts]			# Points rejected
+int	npts				# Number of input points
+int	nvars				# Number of input variables
+double	low_reject, high_reject		# Rejection thresholds
+double	grow				# Rejection radius
+int	nreject				# Number of points rejected (output)
+int	newreject			# Number of new points rej. (output)
+
+int	i, j, i_min, i_max, ilast
+double	sigma, low_cut, high_cut, residual
+pointer	sp, residuals
+
+begin
+#	# Debug.
+#	call eprintf (
+#	    "in_deviant: nl=%d, npts=%d, nvars=%d, low=%g, high=%g, grow=%g\n")
+#	    call pargi (nl)
+#	    call pargi (npts)
+#	    call pargi (nvars)
+#	    call parg$t (low_reject)
+#	    call parg$t (high_reject)
+#	    call parg$t (grow)
+
+	# Initialize.
+	nreject = 0
+	newreject = 0
+
+	# If low_reject and high_reject are zero then just return.
+	if ((low_reject == double (0.0)) && (high_reject == double (0.0)))
+	    return
+
+	# Allocate memory for the residuals.
+	call smark (sp)
+	call salloc (residuals, npts, TY_DOUBLE)
+
+	# Compute the residuals.
+	call nlvectord (nl, x, Memd[residuals], npts, nvars)
+	call asubd (y, Memd[residuals], Memd[residuals], npts)
+
+	# Compute the sigma of the residuals.
+	j = 0
+	sigma = double (0.0)
+	do i = 1, npts {
+	    if ((w[i] != double (0.0)) && (rejpts[i] == NO)) {
+		sigma = sigma + Memd[residuals+i-1] ** 2
+		j = j + 1
+	    } else if (rejpts[i] == YES)
+		nreject = nreject + 1
+	}
+
+	# If there are less than five points for the sigma calculation,
+	# just return.
+
+	if (j < 5) {
+	    call sfree (sp)
+	    return
+	} else
+	    sigma = sqrt (sigma / j)
+
+	# Set the lower and upper cut limits according the the sigma value.
+
+	if (low_reject > double (0.0))
+	    low_cut = -low_reject * sigma
+	else
+	    low_cut = -MAX_REAL
+	if (high_reject > double (0.0))
+	    high_cut = high_reject * sigma
+	else
+	    high_cut = MAX_REAL
+
+	# Reject the residuals exceeding the rejection limits.
+	# A for loop is used instead of do because with region
+	# growing we want to modify the loop index.
+
+	for (i = 1; i <= npts; i = i + 1) {
+
+	    # Do not process points with zero weigth or already rejected.
+	    if ((w[i] == double (0.0)) || (rejpts[i] == YES))
+		next
+
+	    # Reject point, and all other points closer than the growing
+	    # factor.
+
+	    residual = Memd[residuals + i - 1]
+	    if ((residual > high_cut) || (residual < low_cut)) {
+
+		# Determine region to reject.
+		i_min = max (1, int (i - grow))
+		i_max = min (npts, int (i + grow))
+
+		# Reject points from the fit and flag them.
+		do j = i_min, i_max {
+		    if ((abs (x[i] - x[j]) <= grow) && (w[j] != double (0.0)) &&
+		        (rejpts[j] == NO)) {
+			rejpts[j] = YES
+		        newreject = newreject + 1
+			ilast = j
+		    }
+		}
+		i = ilast
+	    }
+	}
+
+	# Free memory.
+	call sfree (sp)
+end