include <foc.h>

# rstat -- print statistics about a pair of reseau entries
# This task computes statistics about reseau entries.  If two reseau entries
# are specified, the mean & std dev etc are computed for the difference
# between the two entries.  If only one entry is specified, the statistics
# are computed for that entry alone.
#
# Phil Hodge, 20-May-1990  Task created

procedure rstat()

char	inres1[SZ_FNAME]	# name of reseau table for first entry
char	rentry[SZ_FNAME]	# name of first entry
char	inres2[SZ_FNAME]	# name of reseau table for second entry
char	dentry[SZ_FNAME]	# name of second entry
#--
pointer sp
pointer diffx, diffy		# scratch for difference between reseau pos'n
pointer rp1, rp2		# pointers to reseau tables
pointer r_respt, d_respt	# reseau entry pointer, first & second
real	avg, stddev		# average & standard deviation of x or y
real	median, mad		# median & mean of absolute values of deviations
int	nres			# total number of reseau marks in entry
int	ngood			# number of reseau not including indef pos
int	reseau			# loop index for reseau number
bool	one_table		# both entries from same table?
bool	one_entry		# only one reseau entry specified?
int	rs_readr()

begin
	call smark (sp)

	call clgstr ("inres1", inres1, SZ_FNAME)
	call clgstr ("entry1", rentry, SZ_FNAME)
	call clgstr ("inres2", inres2, SZ_FNAME)
	call clgstr ("entry2", dentry, SZ_FNAME)

	# Open reseau tables.
	call rstat_open (inres1, inres2, rp1, rp2, one_table)

	# Allocate space for the first entry and read it.
	call rs_alloc (rp1, TY_REAL, r_respt)
	if (rs_readr (rp1, r_respt, rentry, 0) == RES_F_NONEXTENTRY)
	    call error (1, "entry1 not found")

	# Now the second entry.
	if (dentry[1] == EOS || dentry[1] == ' ') {

	    one_entry = true

	} else {

	    call rs_alloc (rp2, TY_REAL, d_respt)
	    if (rs_readr (rp2, d_respt, dentry, 0) == RES_F_NONEXTENTRY) {
		call eprintf ("caution:  second reseau entry not found\n")
		one_entry = true

	    } else {
		one_entry = false

		if (RES_NCOLS(rp1) != RES_NCOLS(rp2) ||
		    RES_NROWS(rp1) != RES_NROWS(rp2))
		call error (1, "size of reseau entries must be the same")
	    }
	}

	# Get the differences.
	nres = RES_NCOLS(rp1) * RES_NROWS(rp1)
	call salloc (diffx, nres, TY_REAL)
	call salloc (diffy, nres, TY_REAL)
	ngood = 0				# initial value
	if (one_entry) {
	    # Fill diffx & diffy with x & y coordinates of reseau marks.
	    do reseau = 0, nres-1 {
		if ( ! Memb[RES_FPT(r_respt)+reseau] ) {
		    Memr[diffx+ngood] = Memr[RES_XPT(r_respt)+reseau]
		    Memr[diffy+ngood] = Memr[RES_YPT(r_respt)+reseau]
		    ngood = ngood + 1
		}
	    }
	} else {
	    # Fill diffx & diffy with differences between reseau coordinates.
	    do reseau = 0, nres-1 {
		if ( ! (Memb[RES_FPT(r_respt)+reseau] ||
			Memb[RES_FPT(d_respt)+reseau]) ) {
		    Memr[diffx+ngood] = Memr[RES_XPT(r_respt)+reseau] -
		    			Memr[RES_XPT(d_respt)+reseau]
		    Memr[diffy+ngood] = Memr[RES_YPT(r_respt)+reseau] -
					Memr[RES_YPT(d_respt)+reseau]
		    ngood = ngood + 1
		}
	    }
	}

	# Get & print the statistics for x & y.

	if (one_table) {
	    call printf ("# table `%s', ")
		call pargstr (inres1)
	} else {
	    call printf ("# tables `%s', `%s', ")
		call pargstr (inres1)
		call pargstr (inres2)
	}
	if (one_entry) {
	    call printf ("entry `%s'\n")
		call pargstr (rentry)
	} else {
	    call printf ("entries `%s', `%s'\n")
		call pargstr (rentry)
		call pargstr (dentry)
	}

	call printf (
"#%5d            mean          stddev          median             mad\n")
	    call pargi (ngood)

	call rstat_avg (Memr[diffx], ngood, avg, stddev)
	call rstat_med (Memr[diffx], ngood, median, mad)
	call printf ("    x: %15.7g %15.7g %15.7g %15.7g\n")
	    call pargr (avg)
	    call pargr (stddev)
	    call pargr (median)
	    call pargr (mad)

	call rstat_avg (Memr[diffy], ngood, avg, stddev)
	call rstat_med (Memr[diffy], ngood, median, mad)
	call printf ("    y: %15.7g %15.7g %15.7g %15.7g\n")
	    call pargr (avg)
	    call pargr (stddev)
	    call pargr (median)
	    call pargr (mad)

	if ( ! one_table )
	    call rs_close (rp2)

	call rs_close (rp1)
	call sfree (sp)
end

# rstat_open -- open reseau tables
# This routine opens the two input reseau tables.

procedure rstat_open (inres1, inres2, rp1, rp2, one_table)

char	inres1[ARB]	# i: name of reseau table for first entry
char	inres2[ARB]	# io: second reseau table (leading blanks removed)
pointer rp1		# o: pointer to struct for inres1
pointer rp2		# o: pointer to struct for inres2; this will be the
			# same as rp1 if one_table is true
bool	one_table	# o: true if inres2 is same as inres1
#--
pointer rs_open()
bool	streq()

begin
	# Remove leading blanks in case inres2 is blank when the user
	# intended it to be null.
	call xt_stripwhite (inres2)

	if (inres2[1] == EOS || streq (inres1, inres2)) {
	    # Get both from the same input table.
	    one_table = true
	    rp1 = rs_open (inres1, READ_ONLY, NULL)
	    rp2 = rp1
	} else {
	    # Separate tables.
	    one_table = false
	    rp1 = rs_open (inres1, READ_ONLY, NULL)
	    rp2 = rs_open (inres2, READ_ONLY, NULL)
	}
end

procedure rstat_avg (x, n, avg, stddev)

real	x[ARB]		# i: array to be averaged
int	n		# i: size of array
real	avg		# o: average of x
real	stddev		# o: standard deviation of x
#--
double	sum, diff
int	k

begin
	if (n < 1) {
	    avg = INDEF
	    stddev = INDEF

	} else if (n == 1) {
	    avg = x[1]
	    stddev = INDEF

	} else {
	    sum = 0.d0

	    do k = 1, n
		sum = sum + x[k]

	    avg = sum / n

	    sum = 0.d0
	    do k = 1, n {
		diff = (x[k] - avg)
		sum = sum + diff * diff
	    }

	    if (n < 2)
		stddev = INDEF
	    else
		stddev = sqrt (sum / (n-1))
	}
end

# rstat_med -- get the median of an array of reals

procedure rstat_med (x, npts, median, mad)

real	x[npts]		# i: input array
int	npts		# i: number of points
real	median		# o: median of a
real	mad		# o: mean of the absolute values of the deviations
#--
pointer	sp, index, ptr
real	sum
int	i, j

int	rstat_compare()
extern	rstat_compare

begin
	if (npts == 1) {
	    median = x[1]
	    mad = 0.
	    return
	} else if (npts == 2) {
	    median = (x[1] + x[2]) / 2
	    mad = abs (x[1] - x[2]) / 2
	    return
	}

	call smark (sp)
	call salloc (index, npts, TY_INT)
	call salloc (ptr, npts, TY_REAL)

	do i = 1, npts
	    Memi[index+i-1] = ptr + i - 1

	call amovr (x, Memr[ptr], npts)

	call qsort (Memi[index], npts, rstat_compare)

	# If npts is odd, i and j will be equal; otherwise, average
	# the middle two points.
	i = (npts + 1) / 2 - 1		# zero indexed
	j = (npts + 2) / 2 - 1
	median = (Memr[Memi[index+i]] + Memr[Memi[index+j]]) / 2

	# Get mad.
	sum = 0.
	do i = 1, npts
	    sum = sum + abs (x[i] - median)
	mad = sum / npts

	call sfree (sp)
end

# compare to give increasing order
int	procedure rstat_compare (i, j)

pointer	i, j		# Array indices to be compared.

begin
	if (Memr[i] < Memr[j])
	    return (-1)
	else if (Memr[i] > Memr[j])
	    return (1)
	else
	    return (0)
end