include <foc.h>

define	RG_AVERAGE	1	# just average the values
define	RG_CLIP		2	# average, clipping outliers
define	RG_MEDIAN	3	# take the median

# resavg -- average entries in a reseau table
# This task takes the median of the reseau positions in a list of entries
# in a reseau table and writes the result as a new entry, either in the
# input table or in an output table.  The output table may be new or may
# already exist.
#
# Phil Hodge, 1-Dec-1989  Task created
# Dave Baxter, Dec-1989  Averaging subroutines written and inserted
# Dave Baxter, 7-Mar-1990 Amended so that reseaux with no measured positions
#                         are correctly dealt with
# Phil Hodge, 20-Feb-1991  Change declaration of streq from int to bool.
# Phil Hodge, 30-Mar-1992  Change letter for "tracking" from I to A.
# Phil Hodge, 15-Jul-1993  In ravg_open, remove ",0,0" from call to tbtacc.

procedure resavg()

char	inres[SZ_FNAME]		# name of input reseau table
char	outres[SZ_FNAME]	# name of output reseau table
char	ientry[SZ_FNAME]	# template entry name in input table
char	oentry[SZ_FNAME]	# name of entry in output table
char	avg_option[SZ_FNAME]	# what kind of average to take
int	option			# what kind of average to take
bool    setindef		# set `bad' reseau positions to INDEF
real    devlim			# mean scatter allowed for valid postions
bool	verbose			# print results?
#--
pointer rlp			# scratch for array of entry pointers
pointer xp			# scr for small array of x coordinates
pointer yp			# scr for small array of y coordinates
pointer irp, orp		# pointers to input, output reseau tables
pointer orespt			# output reseau entry pointer

real	xavg, yavg		# average x & y of a reseau mark
real	rmad			# mean average deviations
real    total			# total mean radial deviations
real    mrmad			# mean deviation for all reseaux
int     nindef			# number of positions set to INDEF
int     counter 		# number of valid reseau positions

int	rlist			# entry template pointer
int	nrows			# number of rows in input table
int	num_ent			# number of entries given by template ientry
int	n			# number of non-null entries to be averaged
int	reseau			# loop index for reseau number
int	nres			# total number of reseau marks in an entry
bool	inplace			# output reseau table same as input?
bool    valid			# bad reseau flag
pointer rs_template()
int	rs_num_rows()
int	clgwrd()
real    clgetr()
bool	clgetb()
string	dictionary	"|average|clip|median|"

begin
	call clgstr ("inres", inres, SZ_FNAME)
	call clgstr ("outres", outres, SZ_FNAME)
	call clgstr ("ientry", ientry, SZ_FNAME)
	call clgstr ("oentry", oentry, SZ_FNAME)
	option = clgwrd ("option", avg_option, SZ_FNAME, dictionary)
	setindef = clgetb ("setindef")
	devlim = clgetr ("devlim")
	verbose = clgetb ("verbose")


	# Open input and output reseau tables.
	call ravg_open (inres, outres, irp, orp, inplace)

	# Open template for input reseau entries.
	rlist = rs_template (irp, ientry)

	nrows = rs_num_rows (irp)	# total number of rows
	if (nrows <= 0)
	    call error (1, "no rows in input table")

	# Allocate space for array of pointers to entry structs and space
	# for values to be averaged.  nrows + 1 elements are allocated for
	# reseau entry pointers; the extra one is for the output entry.

	call malloc (rlp, nrows+1, TY_STRUCT)	# array of pointers
	call malloc (xp, nrows, TY_REAL)	# x values to be averaged
	call malloc (yp, nrows, TY_REAL)	# y values to be averaged

	# Allocate space for each input reseau entry.  The number num_ent
	# of entries that match the template may be fewer than nrows.

	call printf ("#-------------------------\n")
	call printf ("#  Input Reseau Table was:                %-s\n")
		call pargstr ( inres )
	call printf ("#  Chosen method was     :                %-s\n")
		call pargstr ( avg_option )
	call printf ("#-------------------------\n")


	call ravg_alloc (irp, rlist, nrows, Memi[rlp], num_ent, orespt)

	if ( setindef )  {
	call printf ("\n#  setindef = `yes'.   Deviation limit =%7.3f\n")
			call pargr (devlim)
		} else 
		call printf ("\n#  setindef = `no'.\n")
 	call printf ("#----------------------------------------------\n")
		
	# Average x & y coordinates of each reseau mark.
	    if ( verbose )  {
        call printf("\n#Res No.    <X>        <Y>        <d>     Number\n\n")
	    }
	nres = RES_NCOLS(irp) * RES_NROWS(irp)

	counter = 0
	total = 0.0

	do reseau = 0, nres-1 {			# note:  zero indexed

	    # copy coordinates to scratch
	    call ravg_copy (Memi[rlp], num_ent, reseau, Memr[xp], Memr[yp], n)

	    if (n > 0) {

		# average x & y coords of current reseau mark

		valid = true

 		call ravg_average (Memr[xp], Memr[yp], n, option, valid,
			devlim, xavg, yavg, rmad)

		if ( ! setindef )
			valid = true

		if ( valid ) {
			total = total + rmad
			counter = counter + 1
 			Memb[RES_FPT(orespt)+reseau] = false
		}
	    } else {

		# current reseau mark is undefined
		
		valid = false

		Memr[RES_XPT(orespt)+reseau] = INDEFR
		Memr[RES_YPT(orespt)+reseau] = INDEFR
		Memb[RES_FPT(orespt)+reseau] = true

		xavg = INDEFR
		yavg = INDEFR
		rmad = 0.0
	    }

	    if (verbose) {
		call printf ("  %3d    %7.2f    %7.2f    %7.2f      %2d")
		    call pargi (reseau+1)
		    call pargr (xavg)
		    call pargr (yavg)
		    call pargr (rmad)
		    call pargi (n)
		if ( valid )
		    call printf ("\n")
		else
		    call printf ("    set to `INDEF'\n")
	    }
	    if ( ! valid ) {
  		xavg = INDEFR
  		yavg = INDEFR
		Memb[RES_FPT(orespt)+reseau] = true
	    }

	    Memr[RES_XPT(orespt)+reseau] = xavg
	    Memr[RES_YPT(orespt)+reseau] = yavg

	}
	mrmad = total/counter
        nindef = nres-counter

	call clputr ("mdevfmt", mrmad)
	call clputi ("nindef", nindef)

	call printf ("\n Mean Scatter for Format = %7.3f.\n")
	    call pargr ( mrmad )
	call printf ("  %3d invalid positions.\n")
	    call pargi ( nindef )

	# Save averages in output reseau table.

	call rs_t_update (orespt, "A", SZ_R_TRACKING)
	call strcpy (oentry, RES_ENT_ENTRY(orespt), SZ_R_ENTRY)
	call rs_write (orp, orespt, 0)

	if ( ! inplace )
	    call rs_close (orp)

	call rs_close (irp)

	call mfree (yp, TY_REAL)
	call mfree (xp, TY_REAL)
	call mfree (rlp, TY_STRUCT)
end

# ra_open -- open reseau tables

# This routine opens the input and output reseau tables.

procedure ravg_open (inres, outres, irp, orp, inplace)

char	inres[ARB]	# i: name of input reseau table
char	outres[ARB]	# io: output reseau table (leading blanks removed)
pointer irp		# o: pointer to struct for inres
pointer orp		# o: pointer to struct for outres; this will be the
			# same as irp if inplace is true
bool	inplace		# o: true if outres is same as inres
#--
pointer rs_open()
int	tbtacc()
bool	streq()

begin
	call xt_stripwhite (outres)		# remove leading blanks

	if (outres[1] == EOS || streq (inres, outres)) {

	    # Write the result into the input reseau table.

	    inplace = true
	    irp = rs_open (inres, READ_WRITE, NULL)
	    orp = irp

	} else {

	    # Write the result into a different output table.

	    inplace = false
	    irp = rs_open (inres, READ_ONLY, NULL)
	    if (tbtacc (outres) == YES)			# does it exist?
		orp = rs_open (outres, READ_WRITE, NULL)
	    else
		orp = rs_open (outres, NEW_COPY, irp)
	}
end

procedure ravg_alloc (irp, rlist, nrows, irespt, num_ent, orespt)

pointer irp		# i: pointer to input reseau table struct
pointer rlist		# i: pointer to template for input entries
int	nrows		# i: max number of entries
pointer irespt[ARB]	# o: pointers to input entry structures
int	num_ent		# o: actual number of elements in irespt
pointer orespt		# o: pointer to output entry struct
#--
int	k		# loop index
bool	done		# loop-termination flag
pointer rs_readr()
int	rs_list()

begin
	k = 1
	done = false
	while ( ! done ) {
	    call rs_alloc (irp, TY_REAL, irespt[k])
	    if (rs_readr (irp, irespt[k], "", rs_list (rlist)) ==
					RES_F_NONEXTENTRY)
		done = true
	    else
		k = k + 1

	    if (k > nrows+1)
		call error (1, "ravg_alloc:  too many entries found")
	}

	# One more entry pointer was allocated than is needed for input;
	# use it for the output entry.

	orespt = irespt[k]
	num_ent = k - 1			# the last one is invalid

	call printf ("\n#  The following entries were processed:\n\n")

	do k = 1, num_ent {

	    call printf ("                    %-s\n")
		call pargstr ( RES_ENT_ENTRY(irespt[k]) )
	}

	# Assign the numbers of rows and columns and the sfield stuff.

	RES_ENT_NROWS(orespt) = RES_ENT_NROWS(irespt[1])
	RES_ENT_NCOLS(orespt) = RES_ENT_NCOLS(irespt[1])
	RES_ENT_SFIELD1(orespt) = RES_ENT_SFIELD1(irespt[1])
	RES_ENT_SFIELD2(orespt) = RES_ENT_SFIELD2(irespt[1])
	RES_ENT_SFIELD3(orespt) = RES_ENT_SFIELD3(irespt[1])
	RES_ENT_SFIELD4(orespt) = RES_ENT_SFIELD4(irespt[1])
end


# ravg_copy -- copy coordinates to scratch
# This routine copies x and y coordinates of a given reseau mark
# to an output array.  Only good values (i.e. "found" reseau) are copied.

procedure ravg_copy (irespt, num_ent, reseau, x, y, n)

pointer irespt[num_ent]		# i: pointers to entry structures
int	num_ent			# i: number of input reseau entries
int	reseau			# i: current reseau number (zero indexed)
real	x[ARB]			# o: array of x coordinate values
real	y[ARB]			# o: array of y coordinate values
int	n			# o: size of arrays x & y
#--
int	k			# loop index

begin
	n = 0			# initial value
	do k = 1, num_ent {
	    if ( ! Memb[RES_FPT(irespt[k])+reseau] ) {
		n = n + 1
		x[n] = Memr[RES_XPT(irespt[k])+reseau]
		y[n] = Memr[RES_YPT(irespt[k])+reseau]
	    }
	}
end


# ravg_average -- average values
# This routine takes the average or median (or whatever) of an array of values.

procedure ravg_average (x, y, n, option, valid, devlim, xavg, yavg, rmad)

real	x[n]		# i: array of x values
real	y[n]		# i: array of y values
int	n		# i: size of array
int	option		# i: specifies what type of average to take
bool    valid           # io: valid reseau flag
real    devlim		# i: max scatter for valid reseaux
real	xavg		# o: average of x values
real	yavg		# o: average of y values
real	rmad		# o: mean of the absolute values of the deviations
#---

real    sum1
real    xdev, ydev, rdev
int     i

begin
	if (option == RG_AVERAGE) {
	    call ravg_avg (x, y, n, xavg, yavg, rmad)

	} else if (option == RG_CLIP)  {
	    call ravg_avg (x, y, n, xavg, yavg, rmad)
	    call ravg_clip (x, y, n, xavg, yavg, rmad)

	} else if (option == RG_MEDIAN)  {
	    call ravg_med (x, n, xavg)
	    call ravg_med (y, n, yavg)

	    sum1 = 0.
	    do i = 1, n {
		xdev = abs (x[i] - xavg)
		ydev = abs (y[i] - yavg)
		rdev = sqrt (xdev*xdev + ydev*ydev)
		sum1 = sum1 + rdev
	    }

	    rmad = sum1 / n
	}

	if ( rmad > devlim )
	    valid = false
end

procedure ravg_avg (x, y, n, xavg, yavg, rmad)

real	x[n]		# i: array of x values
real	y[n]		# i: array of y values
int	n		# i: size of array
real	xavg		# o: average of x values
real	yavg		# o: average of y values
real	rmad		# o: mean of the absolute values of the deviations
#--
real	sum1, sum2
real    xdev, ydev, rdev
int	i

begin
	sum1 = 0.
	sum2 = 0.
	do i = 1, n {
	    sum1 = sum1 + x[i]
	    sum2 = sum2 + y[i]
	}

	
	xavg = sum1 / n
	yavg = sum2 / n

	sum1 = 0.
	do i = 1, n {
	    xdev = abs (x[i] - xavg)
	    ydev = abs (y[i] - yavg)
	    rdev = sqrt (xdev*xdev + ydev*ydev)
	    sum1 = sum1 + rdev
	}

	rmad = sum1 / n
end


procedure ravg_clip (x, y, n, xavg, yavg, rmad)

real	x[n]		# i: array of x values
real	y[n]		# i: array of y values
int	n		# i: size of array
real	xavg		# i/o: average of x values
real	yavg		# i/o: average of y values
real	rmad		# i/o: mean of the absolute values of the deviations
#--
real	sumr, sumx, sumy
real 	xdev, ydev, rdev
real	xavg2, yavg2
int	i, count

begin

	sumx = 0.0
	sumy = 0.0
	count = 0
	do i = 1, n {
	    xdev = abs (x[i] - xavg)
	    ydev = abs (y[i] - yavg)
	    rdev = sqrt (xdev*xdev + ydev*ydev)
	    if ( rdev < 1.5*rmad )  {
		sumx = sumx + x[i]
		sumy = sumy + y[i]
		count = count + 1
	    }
	}
	xavg2 = sumx / count
	yavg2 = sumy / count

	sumr = 0.0
	count = 0
	do i = 1, n {	
	    xdev = abs (x[i] - xavg)
	    ydev = abs (y[i] - yavg)
	    rdev = sqrt (xdev*xdev + ydev*ydev)
	    if ( rdev < 1.5*rmad )  {
        	xdev = abs (x[i] - xavg2)
   	 	ydev = abs (y[i] - yavg2)
		rdev = sqrt (xdev*xdev + ydev*ydev)
		sumr = sumr + rdev
		count = count + 1
	    }
	}
	xavg = xavg2
	yavg = yavg2
	rmad = sumr / count
	n = count
end


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

procedure ravg_med (x, npts, median)

real	x[npts]		# i: input array
int	npts		# i: number of points
real	median		# o: median of a
#--
int	i, j
pointer	sp, index, ptr

int	ravg_compare()
extern	ravg_compare

begin
	if (npts == 1) {
	    median = x[1]
	    return
	} else if (npts == 2) {
	    median = (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, ravg_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

	call sfree (sp)
end

# compare to give increasing order
int	procedure ravg_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