diff options
Diffstat (limited to 'pkg/utilities/nttools/tstat')
| -rw-r--r-- | pkg/utilities/nttools/tstat/mkpkg | 13 | ||||
| -rw-r--r-- | pkg/utilities/nttools/tstat/thistogram.h | 8 | ||||
| -rw-r--r-- | pkg/utilities/nttools/tstat/thistogram.x | 348 | ||||
| -rw-r--r-- | pkg/utilities/nttools/tstat/thoptions.x | 343 | ||||
| -rw-r--r-- | pkg/utilities/nttools/tstat/tstat.x | 465 |
5 files changed, 1177 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/tstat/mkpkg b/pkg/utilities/nttools/tstat/mkpkg new file mode 100644 index 00000000..f708c769 --- /dev/null +++ b/pkg/utilities/nttools/tstat/mkpkg @@ -0,0 +1,13 @@ +# Update the thistogram application code in the ttools package library +# Author: Phil Hodge, 2-DEC-1988 + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + thistogram.x <error.h> <tbset.h> "thistogram.h" + thoptions.x "thistogram.h" + tstat.x <tbset.h> + ; diff --git a/pkg/utilities/nttools/tstat/thistogram.h b/pkg/utilities/nttools/tstat/thistogram.h new file mode 100644 index 00000000..5c2c0910 --- /dev/null +++ b/pkg/utilities/nttools/tstat/thistogram.h @@ -0,0 +1,8 @@ +# thistogram.h +define NPAR 6 # there are six parameters altogether +define NBINS 1 # got[NBINS] = true means we have a value for nbins +define VLOW 2 +define VHIGH 3 +define DX 4 +define CLOW 5 +define CHIGH 6 diff --git a/pkg/utilities/nttools/tstat/thistogram.x b/pkg/utilities/nttools/tstat/thistogram.x new file mode 100644 index 00000000..7470e1a5 --- /dev/null +++ b/pkg/utilities/nttools/tstat/thistogram.x @@ -0,0 +1,348 @@ +include <error.h> # defines EA_WARN +include <fset.h> # to check whether input or output is redirected +include <tbset.h> +include "thistogram.h" # defines NPAR, etc. + +define MAX_RANGES (SZ_LINE/2) # max number of ranges of row numbers + +# thistogram -- make a histogram of a table column +# +# Phil Hodge, 2-Dec-1988 Task created. +# Phil Hodge, 12-Jan-1989 th_mk_hist: ignore values that are out of range +# Phil Hodge, 17-Mar-1994 Include parameters dx, clow, chigh. +# Phil Hodge, 3-Oct-1995 Modify to use tbn instead of fnt. +# Phil Hodge, 8-Apr-1999 Call tbfpri. +# Phil Hodge, 8-Jun-1999 Set input/output to STDIN/STDOUT if redirected. + +procedure thistogram() + +pointer inlist, outlist # scr for input & output lists of names +char colname[SZ_COLNAME] # column name +int t_nbins # number of bins +double t_vlow, t_vhigh # lower & upper limits for histogram +double t_dx # bin width +double t_clow, t_chigh # centers of first and last bins +char outcolx[SZ_COLNAME] # column name for indep var for histogram +char outcoly[SZ_COLNAME] # column name for dependent var for histogram +#-- +pointer sp +pointer itp, otp # ptr to table descriptor +pointer cptr # ptr to column descriptor +pointer ocpx, ocpy # ptr to col descr for output columns +pointer intab, outtab # scr for names of input & output tables +pointer range_string # string which gives ranges of row numbers +pointer val, counts # scr for histogram: indep & dep var + +# These six parameters are copied from t_... in each loop. +int nbins +double vlow, vhigh, dx, clow, chigh + +pointer list1, list2 +int i, junk +int nrows # number of rows included and not INDEF +int phu_copied # set by tbfpri and ignored +bool listout # is the output ASCII rather than a table? +bool got[NPAR] # flags to specify what we have got +bool find_datamin # true if we need to find minimum data value +bool find_datamax # true if we need to find maximum data value +pointer tbtopn() +double clgetd() +int fstati() +pointer tbnopen() +int clgeti(), tbnget(), tbnlen() +bool streq() + +begin + # Allocate scratch for lists of names and for table names. + call smark (sp) + call salloc (inlist, SZ_FNAME, TY_CHAR) + call salloc (outlist, SZ_FNAME, TY_CHAR) + call salloc (intab, SZ_FNAME, TY_CHAR) + call salloc (outtab, SZ_FNAME, TY_CHAR) + call salloc (range_string, SZ_FNAME, TY_CHAR) + + # Get task parameters. + + if (fstati (STDIN, F_REDIR) == YES) + call strcpy ("STDIN", Memc[inlist], SZ_FNAME) + else + call clgstr ("intable", Memc[inlist], SZ_FNAME) + + if (fstati (STDOUT, F_REDIR) == YES) + call strcpy ("STDOUT", Memc[outlist], SZ_FNAME) + else + call clgstr ("outtable", Memc[outlist], SZ_FNAME) + + call clgstr ("column", colname, SZ_COLNAME) + + # Some of these six parameters may be INDEF. The "t_" prefix + # means these are task parameters; they are assigned to variables + # (same name but without the "t_") in the loop over tables, and + # those variables may be modified within the loop. + t_nbins = clgeti ("nbins") + t_vlow = clgetd ("lowval") + t_vhigh = clgetd ("highval") + t_dx = clgetd ("dx") + t_clow = clgetd ("clow") + t_chigh = clgetd ("chigh") + + call clgstr ("rows", Memc[range_string], SZ_FNAME) + + listout = streq (Memc[outlist], "STDOUT") # ASCII output? + if ( ! listout ) { + call clgstr ("outcolx", outcolx, SZ_COLNAME) + call clgstr ("outcoly", outcoly, SZ_COLNAME) + } + + if (!IS_INDEF(t_dx)) + if (t_dx <= 0.d0) + call error (1, "dx must not be less than or equal to zero") + + # These parameters are interdependent, so compute what was not + # specified from those that were, as far as we can. + call th_options (t_nbins, t_vlow, t_vhigh, t_dx, t_clow, t_chigh, + got, find_datamin, find_datamax) + + # Expand the input table list. + list1 = tbnopen (Memc[inlist]) + + if ( ! listout ) { + # Expand the output table list. + list2 = tbnopen (Memc[outlist]) + if (tbnlen (list1) != tbnlen (list2)) { + call tbnclose (list1) + call tbnclose (list2) + call error (1, + "Number of input and output tables not the same") + } + } + + # Do for each input table. + while (tbnget (list1, Memc[intab], SZ_FNAME) != EOF) { + + # These may be modified within this loop by th_limits. + nbins = t_nbins + vlow = t_vlow + vhigh = t_vhigh + dx = t_dx + clow = t_clow + chigh = t_chigh + + itp = tbtopn (Memc[intab], READ_ONLY, NULL) + call tbcfnd (itp, colname, cptr, 1) + if (cptr == NULL) { + call tbtclo (itp) + call eprintf ("column not found in %s\n") + call pargstr (Memc[intab]) + if ( ! listout ) # skip next output table + junk = tbnget (list2, Memc[outtab], SZ_FNAME) + next + } + + # Get lower & upper limits for the histogram. + iferr { + call th_limits (itp, cptr, Memc[range_string], + nbins, vlow, vhigh, dx, clow, chigh, + got, find_datamin, find_datamax) + } then { + call erract (EA_WARN) + call eprintf ("Table `%s' will be skipped.\n") + call pargstr (Memc[intab]) + call tbtclo (itp) + next + } + + # Get scratch space for the histogram. + call malloc (val, nbins, TY_DOUBLE) + call malloc (counts, nbins, TY_INT) + + # Make the histogram. + call th_mk_hist (itp, cptr, nbins, vlow, vhigh, dx, + Memc[range_string], Memd[val], Memi[counts], nrows) + + if ( listout ) { + call printf ("# %d rows %s\n") + call pargi (nrows) + call pargstr (Memc[intab]) + do i = 1, nbins { + call printf ("%15.7g %8d\n") + call pargd (Memd[val+i-1]) + call pargi (Memi[counts+i-1]) + } + } else { + + # Create output table & define columns. + junk = tbnget (list2, Memc[outtab], SZ_FNAME) + call tbfpri (Memc[intab], Memc[outtab], phu_copied) + otp = tbtopn (Memc[outtab], NEW_FILE, NULL) + call tbcdef (otp, ocpx, + outcolx, "", "", TY_DOUBLE, 1, 1) + call tbcdef (otp, ocpy, + outcoly, "histogram", "", TY_INT, 1, 1) + call tbtcre (otp) + + # Put info records in the header. + call tbhadt (otp, "intable", Memc[intab]) + call tbhadt (otp, "colname", colname) + call tbhadi (otp, "nrows", nrows) + + # Write the values into the output table, and close it. + call tbcptd (otp, ocpx, Memd[val], 1, nbins) + call tbcpti (otp, ocpy, Memi[counts], 1, nbins) + call tbtclo (otp) + } + call tbtclo (itp) + + call mfree (counts, TY_INT) + call mfree (val, TY_DOUBLE) + } + call tbnclose (list1) + if ( ! listout ) + call tbnclose (list2) + call sfree (sp) +end + +# th_limits -- get limits for histogram +# This routine determines the lower and upper limits of data values for +# making a histogram. If either of the input values v1 or v2 is not INDEF +# (i.e. if it was specified by the user), then that value is returned as +# vlow or vhigh respectively. If either or both are INDEF the minimum and +# maximum values in the table column are gotten, and the minimum and maximum +# are extended a little to ensure that the endpoints are included in the +# histogram. The range is extended by (max - min) / (nbins - 1) / 2 +# on each end. The parameters nbins, vlow, vhigh, and dx may be updated. + +procedure th_limits (tp, cptr, range_str, + nbins, vlow, vhigh, dx, clow, chigh, + got, find_datamin, find_datamax) + +pointer tp # i: ptr to table descriptor +pointer cptr # i: ptr to column descriptor +char range_str[ARB] # i: range of row numbers +int nbins # io: number of bins +double vlow, vhigh # io: lower and upper limits +double dx # io: bin width +double clow, chigh # i: centers of low and high bins +bool got[NPAR] # i: flags to specify what we have got +bool find_datamin # i: true if we need to find minimum data value +bool find_datamax # i: true if we need to find maximum data value +#-- +double value # an element gotten from the table +double vmin, vmax # min & max values in the column +int nrows # number of rows in table +int row # row number +int ranges[3,MAX_RANGES] # ranges of row numbers +int nvalues # returned by decode_ranges and ignored +int stat # returned by get_next_number +bool done +int decode_ranges(), get_next_number(), tbpsta() + +begin + if (find_datamin || find_datamax) { + + # We must determine either the minimum or maximum or both. + + if (decode_ranges (range_str, ranges, MAX_RANGES, nvalues) != OK) + call error (1, "bad range of row numbers") + nrows = tbpsta (tp, TBL_NROWS) + + # First get initial values for min & max in column. We can't just + # take the first value because it might be INDEF. + row = 0 # initialize get_next_number + stat = get_next_number (ranges, row) # get first row number + done = (stat == EOF) || (row > nrows) + while ( ! done ) { + call tbegtd (tp, cptr, row, value) + if ( IS_INDEFD(value) ) { + # get next row number + stat = get_next_number (ranges, row) + if ((stat == EOF) || (row > nrows)) + call error (1, "all values are INDEF") + } else { + vmin = value + vmax = value + done = true + } + } + + # Update min & max values. + stat = get_next_number (ranges, row) # get next row number + done = (stat == EOF) || (row > nrows) + while ( ! done ) { + call tbegtd (tp, cptr, row, value) + if ( !IS_INDEFD(value) ) { + if (value < vmin) + vmin = value + if (value > vmax) + vmax = value + } + stat = get_next_number (ranges, row) # get next row number + if ((stat == EOF) || (row > nrows)) + done = true + } + + # Update parameter values. + call th_update (vmin, vmax, nbins, vlow, vhigh, dx, clow, chigh, + got, find_datamin, find_datamax) + } +end + +# th_mk_hist -- make the histogram + +procedure th_mk_hist (tp, cptr, nbins, vlow, vhigh, dx, range_str, + val, counts, nrows) + +pointer tp # i: ptr to table descriptor +pointer cptr # i: ptr to column descriptor +int nbins # i: number of bins +double vlow # i: min value for histogram +double vhigh # i: max value for histogram +double dx # i: bin width +char range_str[ARB] # i: range of row numbers +double val[nbins] # o: array of values at center of bins +int counts[nbins] # o: the histogram, array of counts within bins +int nrows # o: number of rows in range_str and value within limits +#-- +double value # an element gotten from the table +int bin # bin number for output +int totalrows # total number of rows in table +int row # row number +int i +int ranges[3,MAX_RANGES] # ranges of row numbers +int nvalues # returned by decode_ranges and ignored +int stat # returned by get_next_number +bool done +int decode_ranges(), get_next_number(), tbpsta() + +begin + totalrows = tbpsta (tp, TBL_NROWS) + + # Initialize the range of row numbers. + if (decode_ranges (range_str, ranges, MAX_RANGES, nvalues) != OK) + call error (1, "bad range of row numbers") + + do i = 1, nbins { + val[i] = vlow + (i - 0.5d0) * dx # value at center of bin + counts[i] = 0 # initialize histogram + } + + nrows = 0 # initialize counter + row = 0 # initialize get_next_number + stat = get_next_number (ranges, row) # get first row number + done = (stat == EOF) || (row > totalrows) + while ( ! done ) { + call tbegtd (tp, cptr, row, value) + if ( ! IS_INDEF(value) ) { + if (value >= vlow && value < vhigh) { + if (dx > 0.d0) + bin = int ((value - vlow) / dx) + 1 + else + bin = 1 + counts[bin] = counts[bin] + 1 + nrows = nrows + 1 + } + } + stat = get_next_number (ranges, row) # get next row number + done = (stat == EOF) || (row > totalrows) + } +end diff --git a/pkg/utilities/nttools/tstat/thoptions.x b/pkg/utilities/nttools/tstat/thoptions.x new file mode 100644 index 00000000..5b9ce639 --- /dev/null +++ b/pkg/utilities/nttools/tstat/thoptions.x @@ -0,0 +1,343 @@ +include "thistogram.h" # defines NPAR, etc. + +# This file contains th_options and th_update. +# +# Phil Hodge, 18-Mar-1994 Subroutines created. + +# th_options -- different options for specifying limits +# There are different ways to specify the limits and bin spacing. This +# routine computes some parameters given others, based on the following: +# +# dx = (vhigh - vlow) / nbins +# dx = (chigh - clow) / (nbins - 1) +# clow = vlow + dx / 2 +# chigh = vhigh - dx / 2 +# +# Note that vlow and vhigh correspond to task parameters lowval and highval. + +procedure th_options (nbins, vlow, vhigh, dx, clow, chigh, + got, find_datamin, find_datamax) + +int nbins # io: number of bins +double vlow, vhigh # io: lower and upper limits +double dx # io: bin width +double clow, chigh # io: centers of low and high bins +bool got[NPAR] # o: flags to specify what we have got +bool find_datamin # o: true if we need to find minimum data value +bool find_datamax # o: true if we need to find maximum data value +#-- +bool fp_equald() + +begin + # These flags will be reset below if we determine the values. + got[NBINS] = !IS_INDEFI(nbins) + got[VLOW] = !IS_INDEFD(vlow) + got[VHIGH] = !IS_INDEFD(vhigh) + got[DX] = !IS_INDEFD(dx) + got[CLOW] = !IS_INDEFD(clow) + got[CHIGH] = !IS_INDEFD(chigh) + + # Check whether low value is greater than high value. + if (got[VLOW] && got[VHIGH]) + if (vlow > vhigh) + call error (1, "lowval must not be larger than highval") + if (got[CLOW] && got[CHIGH]) + if (clow > chigh) + call error (1, "clow must not be larger than chigh") + + # Further checking. + if (got[VLOW] && got[CLOW]) + if (vlow > clow) + call error (1, "lowval must not be larger than clow") + if (got[CHIGH] && got[VHIGH]) + if (chigh > vhigh) + call error (1, "chigh must not be larger than highval") + if (got[VLOW] && got[CHIGH]) + if (vlow > chigh) + call error (1, "lowval must not be larger than chigh") + if (got[CLOW] && got[VHIGH]) + if (clow > vhigh) + call error (1, "clow must not be larger than highval") + + # Set flags to specify what (if anything) we must get from the data. + # These may be reset below. + find_datamin = (!got[VLOW] && !got[CLOW]) + find_datamax = (!got[VHIGH] && !got[CHIGH]) + + if (got[DX]) { + + # Was the lower limit specified by the user? + if (got[CLOW]) { + if (got[VLOW]) { + if (!fp_equald (clow, vlow + dx/2.d0)) + call error (1, "values of dx, clow, lowval conflict") + } else { + vlow = clow - dx / 2.d0 + got[VLOW] = true + } + } else if (got[VLOW]) { + clow = vlow + dx / 2.d0 + got[CLOW] = true + } + + # Was the upper limit specified? + if (got[CHIGH]) { + if (got[VHIGH]) { + if (!fp_equald (chigh, vhigh - dx/2.d0)) + call error (1, "values of dx, chigh, highval conflict") + } else { + vhigh = chigh + dx / 2.d0 + got[VHIGH] = true + } + } else if (got[VHIGH]) { + chigh = vhigh - dx / 2.d0 + got[CHIGH] = true + } + + # Was the number of bins specified? + if (got[NBINS]) { + if (got[VLOW] && got[VHIGH]) { + if (!fp_equald (vhigh - vlow, dx * nbins)) + call error (1, "specified values for limits conflict") + } else if (got[VLOW]) { + vhigh = vlow + dx * nbins + chigh = vhigh - dx / 2.d0 + got[VHIGH] = true + got[CHIGH] = true + find_datamax = false + } else if (got[VHIGH]) { + vlow = vhigh - dx * nbins + clow = vlow + dx / 2.d0 + got[VLOW] = true + got[CLOW] = true + find_datamin = false + } + + } else if (got[VLOW] && got[VHIGH]) { + nbins = nint ((vhigh - vlow) / dx) + if (nbins * dx < vhigh - vlow) + nbins = nbins + 1 # round up + got[NBINS] = true + } + + } else if (got[NBINS]) { # but we don't have dx + + if (nbins == 1) { + if (!got[VLOW] && !got[VHIGH] && (got[CLOW] || got[CHIGH])) { + call eprintf ( + "nbins = 1, clow or chigh was specified, but dx was not.\n") + call error (1, "must specify dx for this case") + } + } + + if (got[VLOW] && got[VHIGH]) { + + dx = (vhigh - vlow) / double(nbins) + got[DX] = true + if (got[CLOW]) { + if (!fp_equald (clow, vlow + dx/2.d0)) + call error (1, "clow conflicts with other parameters") + } else { + clow = vlow + dx / 2.d0 + got[CLOW] = true + } + if (got[CHIGH]) { + if (!fp_equald (chigh, vhigh - dx/2.d0)) + call error (1, "chigh conflicts with other parameters") + } else { + chigh = vhigh - dx / 2.d0 + got[CHIGH] = true + } + + } else if (got[CLOW] && got[CHIGH]) { + + if (nbins == 1) { + if (!fp_equald (clow, chigh)) + call error (1, "nbins = 1, but clow != chigh") + } else { + dx = (chigh - clow) / (double(nbins) - 1.d0) + got[DX] = true + if (got[VLOW]) { + if (!fp_equald (vlow, clow - dx/2.d0)) + call error (1, + "lowval conflicts with other parameters") + } else { + vlow = clow - dx / 2.d0 + got[VLOW] = true + } + if (got[VHIGH]) { + if (!fp_equald (vhigh, chigh + dx/2.d0)) + call error (1, + "highval conflicts with other parameters") + } else { + vhigh = chigh + dx / 2.d0 + got[VHIGH] = true + } + } + + } else if (got[CLOW] && got[VLOW]) { + + dx = (clow - vlow) * 2.d0 + vhigh = vlow + dx * nbins + chigh = vhigh - dx / 2.d0 + got[DX] = true + got[VHIGH] = true + got[CHIGH] = true + find_datamax = false + + } else if (got[CHIGH] && got[VHIGH]) { + + dx = (vhigh - chigh) * 2.d0 + vlow = vhigh - dx * nbins + clow = vlow + dx / 2.d0 + got[DX] = true + got[VLOW] = true + got[CLOW] = true + find_datamin = false + + } else if (got[CLOW] && got[VHIGH]) { + + dx = (vhigh - clow) / (double(nbins) - 0.5d0) + vlow = vhigh - dx * nbins + chigh = vhigh - dx / 2.d0 + got[DX] = true + got[VLOW] = true + got[CHIGH] = true + + } else if (got[VLOW] && got[CHIGH]) { + + dx = (chigh - vlow) / (double(nbins) - 0.5d0) + clow = vlow + dx / 2.d0 + vhigh = vlow + dx * nbins + got[DX] = true + got[CLOW] = true + got[VHIGH] = true + + } + + } else if (got[CLOW] && got[VLOW]) { # but neither dx nor nbins + + dx = (clow - vlow) * 2.d0 + got[DX] = true + + } else if (got[CHIGH] && got[VHIGH]) { # but neither dx nor nbins + + dx = (vhigh - chigh) * 2.d0 + got[DX] = true + + } else { + call error (1, "you must specify either nbins or dx (or both)") + } +end + +# th_update -- update the limits +# We now have the minimum and maximum data values from the table, +# so we can fill in the values that the user did not specify. +# We need nbins, vlow, vhigh, and dx. The values of clow and chigh +# are not modified, even if they are still INDEF; the array of flags +# (got) must not be modified, as these flags are used for different +# tables if there is more than one in the input list. + +procedure th_update (vmin, vmax, nbins, vlow, vhigh, dx, clow, chigh, + got, find_datamin, find_datamax) + +double vmin, vmax # i: min and max data values for current table +int nbins # io: number of bins +double vlow, vhigh # io: lower and upper limits +double dx # io: bin width +double clow, chigh # i: centers of low and high bins +bool got[NPAR] # i: flags that specify what parameters we already have +bool find_datamin # i: true if we had to find minimum data value +bool find_datamax # i: true if we had to find maximum data value +#-- +double delta # for expanding the range to include endpoints + +begin + if (find_datamin && find_datamax) { + + # Center the data within the range. We do have either + # dx or nbins or both. + if (got[DX] && got[NBINS]) { + delta = (dx * nbins - (vmax - vmin)) / 2.d0 + } else if (got[NBINS]) { + if (nbins == 1) + delta = (vmax - vmin) / 100.d0 # 100 is arbitrary + else + delta = (vmax - vmin) / (nbins - 1.d0) / 2.d0 + } else if (got[DX]) { + nbins = nint ((vmax - vmin) / dx) + if (nbins * dx <= vmax - vmin) + nbins = nbins + 1 # round up + delta = (dx * nbins - (vmax - vmin)) / 2.d0 + } + vlow = vmin - delta + vhigh = vmax + delta + if (!got[DX]) + dx = (vhigh - vlow) / nbins + + } else if (find_datamin) { + + # We don't have both dx and nbins. If we did, we would have + # calculated vlow from vhigh, dx, and nbins. + if (got[DX]) { + + nbins = nint ((vhigh - vmin) / dx) # vhigh is known + if (nbins * dx <= vhigh - vmin) + nbins = nbins + 1 # round up + vlow = vhigh - dx * nbins + + } else if (got[NBINS]) { + + if (got[VHIGH]) { + if (nbins == 1) + delta = (vhigh - vmin) / 100.d0 # 100 is arbitrary + else + delta = (vhigh - vmin) / (nbins - 1.d0) / 2.d0 + vlow = vmin - delta + dx = (vhigh - vlow) / nbins + } else { # we have chigh but not vhigh + if (nbins == 1) { + vlow = vmin # this case doesn't make sense + vhigh = chigh + (chigh - vmin) + dx = vhigh - vlow + } else { # set clow = vmin + dx = (chigh - vmin) / (nbins - 1.d0) + vlow = vmin - dx / 2.d0 + vhigh = chigh + dx / 2.d0 + } + } + } + + } else if (find_datamax) { + + # For this case as well, we don't have both dx and nbins. + if (got[DX]) { + + nbins = nint ((vmax - vlow) / dx) # vlow is known + if (nbins * dx <= vmax - vlow) + nbins = nbins + 1 # round up + vhigh = vlow + dx * nbins + + } else if (got[NBINS]) { + + if (got[VLOW]) { + if (nbins == 1) + delta = (vmax - vlow) / 100.d0 # 100 is arbitrary + else + delta = (vmax - vlow) / (nbins - 1.d0) / 2.d0 + vhigh = vmax + delta + dx = (vhigh - vlow) / nbins + } else { # we have clow but not vlow + if (nbins == 1) { + vhigh = vmax # this case doesn't make sense + vlow = clow - (vmax - clow) + dx = vhigh - vlow + } else { # set chigh = vmax + dx = (vmax - clow) / (nbins - 1.d0) + vlow = vmin - dx / 2.d0 + vhigh = chigh + dx / 2.d0 + } + } + } + } +end diff --git a/pkg/utilities/nttools/tstat/tstat.x b/pkg/utilities/nttools/tstat/tstat.x new file mode 100644 index 00000000..5ef2b435 --- /dev/null +++ b/pkg/utilities/nttools/tstat/tstat.x @@ -0,0 +1,465 @@ +include <error.h> +include <fset.h> # to check whether input or output is redirected +include <tbset.h> + +define MAX_RANGES (SZ_LINE/2) # max number of ranges of row numbers +define NUM_COL 8 # number of output columns + +# tstat -- get statistics for a table column +# This task gets the mean, standard deviation, median, and minimum & +# maximum values for a table column. +# +# Phil Hodge, 8-Dec-1988 Task created. +# Phil Hodge, 14-Mar-1989 Also compute the median; fix bug in std dev. +# Phil Hodge, 27-May-1992 Print INDEF if no values in range. +# Phil Hodge, 31-Jul-1992 Print column name after table name. +# Phil Hodge, 11-Jan-1992 Use asrtd instead of asokd for median if nr is even. +# Phil Hodge, 3-Oct-1995 Modify to use tbn instead of fnt. +# Phil Hodge, 25-Apr-1997 Use asrtd instead of asokd regardless of nr. +# Phil Hodge, 26-Mar-1998 Get all elements of array columns; ignore the +# column parameter if an input table has only one column. +# Phil Hodge, 8-Jun-1999 Set input/output to STDIN/STDOUT if redirected. +# Phil Hodge, 2-Jan-2001 If the input was redirected, and one command-line +# argument was specified, take that argument to be the +# column name rather than the table name. + +procedure tstat() + +pointer inlist # scratch for list of input table names +char cl_colname[SZ_COLNAME] # column name gotten from cl parameter +double lowlim # lower & upper limits for histogram +double highlim # lower & upper limits for histogram +pointer range_string # string which gives ranges of row numbers +char nam_table[SZ_FNAME] # column name for table name +char nam_name[SZ_COLNAME] # column name for column name +char nam_mean[SZ_COLNAME] # column name for mean +char nam_stddev[SZ_COLNAME] # column name for standard deviation +char nam_med[SZ_COLNAME] # column name for median +char nam_min[SZ_COLNAME] # column name for minimum +char nam_max[SZ_COLNAME] # column name for maximum +char nam_nrows[SZ_COLNAME] # column name for number of good rows +#-- +pointer list1 # for list of input tables +pointer sp +pointer itp, otp # ptr to table descriptor +pointer cptr # ptr to column descriptor +pointer ocp[NUM_COL] # ptrs to col descriptors for output columns +pointer intab, outtab # scr for names of input & output tables +char colname[SZ_COLNAME] # column name +double vmean # mean value +double vstddev # standard deviation of values +double vmedian # median value +double vmin, vmax # minimum & maximum values +int row # output row number +int nrows # number of rows included and not INDEF +int ntables # number of tables in input list +bool listout # ASCII output? +bool tabout # table output? +bool new_table # true if output table does not already exist +int nargs # number of command-line arguments +bool in_redir # is input redirected? +int clgeti() +double clgetd() +int fstati() +pointer tbnopen() +int tbnget(), tbnlen() +pointer tbtopn() +pointer tbcnum() +int tbtacc(), tbpsta() +bool streq() + +begin + # Allocate scratch for lists of names and for table names. + call smark (sp) + call salloc (inlist, SZ_LINE, TY_CHAR) + call salloc (intab, SZ_LINE, TY_CHAR) + call salloc (outtab, SZ_LINE, TY_CHAR) + call salloc (range_string, SZ_LINE, TY_CHAR) + + # Get task parameters. + + nargs = clgeti ("$nargs") + in_redir = fstati (STDIN, F_REDIR) == YES + + if (in_redir) + call strcpy ("STDIN", Memc[inlist], SZ_LINE) + else + call clgstr ("intable", Memc[inlist], SZ_LINE) + + if (fstati (STDOUT, F_REDIR) == YES) + call strcpy ("STDOUT", Memc[outtab], SZ_LINE) + else + call clgstr ("outtable", Memc[outtab], SZ_LINE) + + cl_colname[1] = EOS # initial value + lowlim = clgetd ("lowlim") # these limits may be INDEF + highlim = clgetd ("highlim") + call clgstr ("rows", Memc[range_string], SZ_LINE) + + # ASCII output? table output? create a new output table? + listout = streq (Memc[outtab], "STDOUT") + if ( listout || (Memc[outtab] == EOS) ) { + tabout = false + } else { + tabout = true + new_table = (tbtacc (Memc[outtab]) == NO) + } + + if (tabout) { + call clgstr ("n_tab", nam_table, SZ_FNAME) + call clgstr ("n_nam", nam_name, SZ_COLNAME) + call clgstr ("n_nrows", nam_nrows, SZ_COLNAME) + call clgstr ("n_mean", nam_mean, SZ_COLNAME) + call clgstr ("n_stddev", nam_stddev, SZ_COLNAME) + call clgstr ("n_median", nam_med, SZ_COLNAME) + call clgstr ("n_min", nam_min, SZ_COLNAME) + call clgstr ("n_max", nam_max, SZ_COLNAME) + } + + # Expand the input table list. + list1 = tbnopen (Memc[inlist]) + + ntables = tbnlen (list1) + + if (listout) { + row = 0 # just to have a definite value + if (ntables > 1) { + call printf ("# nrows") + call printf (" mean stddev median") + call printf (" min max\n\n") + } + } else if (tabout) { + # Create output table (or open existing table) & define columns. + if (new_table) + otp = tbtopn (Memc[outtab], NEW_FILE, NULL) + else + otp = tbtopn (Memc[outtab], READ_WRITE, NULL) + + call tbcfnd (otp, nam_table, ocp[1], 1) + call tbcfnd (otp, nam_name, ocp[2], 1) + call tbcfnd (otp, nam_nrows, ocp[3], 1) + call tbcfnd (otp, nam_mean, ocp[4], 1) + call tbcfnd (otp, nam_stddev, ocp[5], 1) + call tbcfnd (otp, nam_med, ocp[6], 1) + call tbcfnd (otp, nam_min, ocp[7], 1) + call tbcfnd (otp, nam_max, ocp[8], 1) + + if (ocp[1] == NULL) + call tbcdef (otp, ocp[1], + nam_table, "", "", -SZ_FNAME, 1, 1) + if (ocp[2] == NULL) + call tbcdef (otp, ocp[2], + nam_name, "", "", -SZ_COLNAME, 1, 1) + if (ocp[3] == NULL) + call tbcdef (otp, ocp[3], + nam_nrows, "", "", TY_INT, 1, 1) + if (ocp[4] == NULL) + call tbcdef (otp, ocp[4], + nam_mean, "", "", TY_DOUBLE, 1, 1) + if (ocp[5] == NULL) + call tbcdef (otp, ocp[5], + nam_stddev, "", "", TY_DOUBLE, 1, 1) + if (ocp[6] == NULL) + call tbcdef (otp, ocp[6], + nam_med, "", "", TY_DOUBLE, 1, 1) + if (ocp[7] == NULL) + call tbcdef (otp, ocp[7], + nam_min, "", "", TY_DOUBLE, 1, 1) + if (ocp[8] == NULL) + call tbcdef (otp, ocp[8], + nam_max, "", "", TY_DOUBLE, 1, 1) + if (new_table) + call tbtcre (otp) + + row = tbpsta (otp, TBL_NROWS) # append to whatever is there + } + + # Do for each input table. + while (tbnget (list1, Memc[intab], SZ_LINE) != EOF) { + + itp = tbtopn (Memc[intab], READ_ONLY, NULL) + if (tbpsta (itp, TBL_NCOLS) == 1) { + cptr = tbcnum (itp, 1) + call tbcigt (cptr, TBL_COL_NAME, colname, SZ_COLNAME) + } else { + if (cl_colname[1] == EOS) { + if (nargs == 1 && in_redir) { + # in this case, the one argument should be column name + call clgstr ("intable", cl_colname, SZ_COLNAME) + # update par file + call clpstr ("intable", "STDIN") + call clpstr ("column", cl_colname) + } else { + call clgstr ("column", cl_colname, SZ_COLNAME) + } + } + call strcpy (cl_colname, colname, SZ_COLNAME) + call tbcfnd (itp, colname, cptr, 1) + } + if (cptr == NULL) { + call tbtclo (itp) + call eprintf ("column %s not found in %s\n") + call pargstr (colname) + call pargstr (Memc[intab]) + next + } + row = row + 1 + + if (listout) { + call printf ("# %s %s\n") + call pargstr (Memc[intab]) + call pargstr (colname) + if (ntables == 1) { + call printf ("# nrows") + call printf (" mean stddev median") + call printf (" min max\n") + } + } + + # Get statistics. + iferr { + call ts_values (itp, cptr, Memc[range_string], lowlim, highlim, + vmean, vstddev, vmedian, vmin, vmax, nrows) + } then { + call erract (EA_WARN) + next + } + + if (listout) { + + call printf ("%5d %17.10g %13.6g %13.6g %13.6g %13.6g\n") + call pargi (nrows) + call pargd (vmean) + call pargd (vstddev) + call pargd (vmedian) + call pargd (vmin) + call pargd (vmax) + + } else if (tabout) { + + # Write the values into the output table. + call tbeptt (otp, ocp[1], row, Memc[intab]) + call tbeptt (otp, ocp[2], row, colname) + call tbepti (otp, ocp[3], row, nrows) + call tbeptd (otp, ocp[4], row, vmean) + call tbeptd (otp, ocp[5], row, vstddev) + call tbeptd (otp, ocp[6], row, vmedian) + call tbeptd (otp, ocp[7], row, vmin) + call tbeptd (otp, ocp[8], row, vmax) + } + + call tbtclo (itp) # close current input table + } + + # Close the output table. It may be that nothing was written to it. + if (tabout) + call tbtclo (otp) + + # Save the results (from the last input table) in cl parameters. + call clputi ("nrows", nrows) + call clputd ("mean", vmean) + call clputd ("stddev", vstddev) + call clputd ("median", vmedian) + call clputd ("vmin", vmin) + call clputd ("vmax", vmax) + + call tbnclose (list1) + call sfree (sp) +end + +# ts_values -- get statistics for a table column +# This routine gets the mean, standard deviation, minimum and maximum +# values for a table column. If lower and/or upper cutoff limits were +# specified (i.e. are not INDEF) then values outside that range will not +# be included in the statistics. INDEF values are also not included. + +procedure ts_values (tp, cptr, range_str, lowlim, highlim, + vmean, vstddev, vmedian, vmin, vmax, nrows) + +pointer tp # i: ptr to table descriptor +pointer cptr # i: ptr to column descriptor +char range_str[ARB] # i: range of row numbers +double lowlim # i: lower cutoff of values to be included +double highlim # i: upper cutoff of values to be included +double vmean # o: mean value +double vstddev # o: standard deviation of values +double vmedian # o: median value +double vmin # o: minimum value +double vmax # o: maximum value +int nrows # o: number of rows included in the statistics +#-- +pointer sp +pointer val # scratch for array of values (for median) +pointer descrip # column selector descriptor +double value # an element gotten from the table +double sum, sumsq # for accumulating sums +double smin, smax # temp min & max +double diff # current value minus first good value +int all_elem # total number of elements in column +int selrows # number of selected rows +int nelem # number of elements in one cell +int nret # number returned (ignored) +int nr # current number of rows +int row # row number +int ranges[3,MAX_RANGES] # ranges of row numbers +int nvalues # returned by decode_ranges and ignored +int stat # returned by get_next_number +int i, j # loop indexes +bool chklow, chkhi # were low (high) limits specified? +bool val_ok # is current value within limits? +bool done # loop-termination flag +int decode_ranges(), get_next_number() +int tbagtd() +errchk tbagtd, tbegtd, tcs_rdaryd + +begin + if (decode_ranges (range_str, ranges, MAX_RANGES, nvalues) != OK) { + call eprintf ("rows = `%s'\n") + call pargstr (range_str) + call error (1, "Range of row numbers is invalid.") + } + + # Get the number of elements per table cell and the number of + # rows selected using row-selector syntax (as opposed to using + # the 'rows' task parameter). + call tbcnel1 (tp, cptr, descrip, nelem, selrows) + + # Find out how many rows there are in the table, restricted by + # either a row selector or the 'rows' task parameter, or both. + if (range_str[1] == '-' || range_str[1] == EOS) { + all_elem = selrows * nelem + } else { + # Count the number of rows specified. + i = 0 # count of row numbers + row = 0 # initialize get_next_number + done = false + while (!done) { + stat = get_next_number (ranges, row) + if (stat == EOF || row > selrows) + done = true + else + i = i + 1 + } + all_elem = i * nelem + } + + # Allocate scratch space to hold an entire column. + call smark (sp) + call salloc (val, all_elem, TY_DOUBLE) + + row = 0 # reinitialize get_next_number + stat = get_next_number (ranges, row) + done = (stat == EOF || row > selrows) + + # Get the data. + i = 1 + while (!done) { + + if (descrip == NULL) { + if (nelem == 1) + call tbegtd (tp, cptr, row, Memd[val+i-1]) + else + nret = tbagtd (tp, cptr, row, Memd[val+i-1], 1, nelem) + } else { + call tcs_rdaryd (tp, descrip, row, all_elem-i+1, + nret, Memd[val+i-1]) + } + + i = i + nelem + + stat = get_next_number (ranges, row) + done = (stat == EOF || row > selrows) + } + if (all_elem != i - 1) + call error (1, "not all elements read from column") + + chklow = ! IS_INDEFD(lowlim) + chkhi = ! IS_INDEFD(highlim) + + # Check which values to include. + sum = 0.d0 + nr = 0 + do i = 0, all_elem-1 { # zero indexed + value = Memd[val+i] + if (!IS_INDEFD(value)) { + val_ok = true # an initial value + if (chkhi && (value > highlim)) { + val_ok = false + Memd[val+i] = INDEFD + } + if (chklow && (value < lowlim)) { + val_ok = false + Memd[val+i] = INDEFD + } + if (val_ok) { + nr = nr + 1 + sum = sum + value + } + } + } + + if (nr < 1) { + # No rows with valid data, so set the output to INDEF. + nrows = 0 + vmean = INDEFD + vstddev = INDEFD + vmedian = INDEFD + vmin = INDEFD + vmax = INDEFD + return + } else { + # Assign two of the output values. + nrows = nr + vmean = sum / double(nr) + } + + # Move good data to the beginning of the array. + if (nr < all_elem) { + j = 0 + do i = 0, all_elem-1 { # zero indexed + if (!IS_INDEFD(Memd[val+i])) { + Memd[val+j] = Memd[val+i] + j = j + 1 + } + } + } + + # Find the min, max, standard deviation of the good values. + sumsq = 0.d0 + smin = Memd[val] + smax = Memd[val] + do i = 0, nr-1 { # zero indexed + value = Memd[val+i] + diff = value - vmean + sumsq = sumsq + diff * diff + if (value < smin) + smin = value + if (value > smax) + smax = value + } + + # Assign output values. + + vmin = smin + vmax = smax + + if (nr > 1) + vstddev = sqrt (sumsq / double(nr-1)) + else + vstddev = INDEFD + + # Determine the median. If there are an even number of values, the + # middle two are averaged. + if (nr < 3) { + vmedian = vmean + } else { + call asrtd (Memd[val], Memd[val], nr) + if (nr / 2 * 2 == nr) { # even number of values? + vmedian = (Memd[val+nr/2-1] + Memd[val+nr/2]) / 2.d0 + } else { + vmedian = Memd[val+nr/2] + } + } + + call sfree (sp) +end |