Table length reduction

OK, I said I'd have a go at a program to intelligently
reduce the size of IBIS tables so that we can have high
resolution where needed and still keep below the 100-point
limit. Here is is, attached. The language will look
very, very strange to most of you; it's Icon, an academic
language that's sort of a cross between C and Prolog.
You can get documentation and a copy for most any platform at
http://www.cs.arizona.edu/icon/

The program as written is pretty crude. It takes JUST the
table data and thins it out using a cubic interpolation.
Command-line syntax is

decimate <pts> <file>

Where <pts> is the maximum number of points in the output,
and <file> is the source. Output is to standard output, and
if no input file is specified it defaults to standard input.
It expects each line to have the same number of fields (and
will barf if they don't match). For unjustifiable reasons it
does NOT accept scientific notation in; it does accept 'munp'
suffixes. The output is scientific. If I'd been more awake
I'd have done both in scientific and left formatting up to a
perl or awk script.

That said, I've experimented it several data sets and it seems
to produce reasonable results. Feedback, improvements, or a
complete rewrite to C welcome.

-- 
D. C. Sessions
dc.sessions@vlsi.com

# Data row: an independent variable and a set of dependent variables
record t_row( x, yvals )

# Get a sequence of words from a line of text,
# where the cset 'delim' separates words
procedure get_word( delim, line )
  local left, right, word
  left := right := 1
  while right <= *line do {
    left := many( delim, line, right ) | right
    right := many( &ascii -- delim, line, left ) | left
    suspend word := line[ left:right ]
    }
  end

# Convert a word to a numerical value.
# The numerical portion consists of the usual
# numerical prefix optionally followed by scale
# of the munp type. Maybe if I get ambitious I'll
# add a scientific notation option.
procedure get_numeric( word )
  local base
  base := get_word( ~( &digits ++ '.-' ), word )
  return case word[ *base+1 ] of {
    "m" : base / 1E3
    "u" : base / 1E6
    "n" : base / 1E9
    "p" : base / 1E12
    default : base
    } | base
  end

# Translate a table into a list of structures
procedure get_fn_set( source )
  local line, lineval, fn
  fn := []
  while line := read( source ) do {
    lineval := []
    every put( lineval, get_numeric( get_word( ' \t', line ) ) )
    put( fn, t_row( lineval[ 1 ], copy( lineval[ 2 : *lineval+1 ] ) ) )
    }
  return fn
  end

# Extract a single function from a list of structures
# In other words, select a single column from the dependent
# variable collection.
procedure get_single_fn( fn_set, column )
  local row, result
  result := []
  every row := !fn_set do
    put( result, t_row( row.x, row.yvals[ column ] ) )
  return result
  end

# Interpolate a value from at 'new_x' based on
# the values of the list of rows 'neighbors'
procedure interpolate( new_x, neighbors )
  local pick, flick, product, result
  result := 0
  every pick := !neighbors do {
    product := pick.yvals
    every ( flick := !neighbors ) ~=== pick do
      product *:= ( new_x - flick.x )/( pick.x - flick.x )
    result +:= product
    }
  return result
  end

# Get an interval of 2n+1 centered at 'ctr' if
# possible but bounded by 1 and 'mx'
procedure left_bound( ctr, n, mx )
  return if ctr < n + 1
    then 1
    else if ctr > mx - n
      then mx - ( n + n )
      else ctr - n
  end

# Select a list of 2* 'n' neighbors for 'row' from list 'fn'
# and return the selected list
procedure pick_window( row, n, fn )
  local lend, rend
  lend := left_bound( row, n, *fn )
  rend := lend + n + n
  return fn[ lend:row ] ||| fn[ row+1:rend+1 ]
  end

# Compare the interpolated value at a row to
# the actual value. Return the relative error
procedure epsilon( row, n, fn, scale )
  local ival
  ival := interpolate( fn[ row ].x, pick_window( row, n, fn ) )
  return abs( ( ival - fn[ row ].yvals )/scale )
  end

# Build a table of error terms from a set of functions
# using an interpolation of order ( 2n-1 )
procedure err_tab( n, fn_set )
  local cn, j, row, result, scale, column
  result := table( 0 )
  every insert( result, !fn_set, 0 )
  every cn := 1 to *(fn_set[1].yvals) do {
    column := get_single_fn( fn_set, cn )
    scale := 0
    every row := !column do scale := scale < abs( row.yvals )
    every j := 1 to *column do
      result[ fn_set[ j ] ] := result[ fn_set[ j ] ] < epsilon( j, 2, column, scale )
    }
  return result
  end

# Decide which table elements to keep and which to cull
# Start by culling the elements which have the least
# error values, but keep their neighbors for this pass
# Return an edited function list
procedure cull( n, fn_set )
  local pick, j, l, r, etab, index_tab, cull_list, cull_tab, fn_size
  fn_size := *fn_set
  etab := err_tab( n, fn_set )
  index_tab := table( 1 )
  cull_tab := []
  every j := 1 to fn_size do {
    insert( index_tab, fn_set[ j ], j )
    put( cull_tab, 0 ) # Start with the decision unknown
    }
  cull_tab[ 1 ] := 1 # keep the first and last elements
  cull_tab[ fn_size ] := 1
                                # sort to a record list by values
  cull_list := sort( etab, 2 )[1:(fn_size/(n+n+1))] # Limit the number of culls
  every pick := index_tab[ (!cull_list)[1] ] do
    if 1 > cull_tab[ pick ] then {
      cull_tab[ pick ] := -1 # Ditch it
      l := left_bound( pick, n, *fn_set )
      r := l + n + n
      every cull_tab[ pick ~= ( l to r ) ] := 1 # Keep the neighbors
      }
  result := []
  every pick := 1 to *fn_set do
    if cull_tab[ pick ] >= 0
      then put( result, fn_set[ pick ] )
# else {
# writes( "Culled " )
# write_fn_line( fn_set[ pick ] )
# }
  return result
  end

# Round off a number 'x' to 'n' digits precision
procedure roundoff( n, x )
  base := log( 10 )
  if x = 0
    then return 0
    else {
      scale := integer( log( abs( x ), 10 ) )
      fudge := integer( exp( ( n - scale )*base ) )
# writes( "(" || n - scale || " " || fudge || ")" )
      return integer( x * fudge ) / real( fudge )
      }
  end

# Write out a line of a function
procedure write_fn_line( fn_line )
  local cwidth
  cwidth := 12
  writes( right( roundoff( 3, fn_line.x ), cwidth, " " ) )
  every writes( right( roundoff( 3, !fn_line.yvals ), cwidth, " " ) )
  write( )
  end

procedure main( argv )
  limit := pop( argv )
  source := open( pop( argv, src ) ) | &input
   fn_set := get_fn_set( source )
   foo := fn_set
   while *foo > limit do foo := cull( 2, foo )
   every bar := !foo do write_fn_line( bar )
  end
Received on Mon May 8 10:32:13 2000