cmd/pzip/pin.c

/***********************************************************************
*                                                                      *
*               This software is part of the ast package               *
*          Copyright (c) 1998-2011 AT&T Intellectual Property          *
*                      and is licensed under the                       *
*                 Eclipse Public License, Version 1.0                  *
*                    by AT&T Intellectual Property                     *
*                                                                      *
*                A copy of the License is available at                 *
*          http://www.eclipse.org/org/documents/epl-v10.html           *
*         (with md5 checksum b35adb5213ca9657e911e9befb180842)         *
*                                                                      *
*              Information and Software Systems Research               *
*                            AT&T Research                             *
*                           Florham Park NJ                            *
*                                                                      *
*                 Glenn Fowler <gsf@research.att.com>                  *
*                                                                      *
***********************************************************************/
#pragma prototyped

/*
 * Glenn Fowler
 * Adam Buchsbaum
 * AT&T Research
 *
 * induce a column partition on fixed row data
 */

static const char usage[] =
"[-?\n@(#)$Id: pin (AT&T Research) 2003-07-17 $\n]"
USAGE_LICENSE
"[+NAME?pin - induce a pzip partition on fixed record data]"
"[+DESCRIPTION?\bpin\b induces a \bpzip\b(1) column partition on data files"
"   of fixed length rows (records) and columns (fields). If a partition"
"   file is specified then that partition is refined. A partition file,"
"   suitable for use by \bpin\b, \bpzip\b(1) or \bpop\b(1) is listed"
"   on the standard output. The input \afile\a is referred to as"
"   \atraining data\a. \b--size\b allows more than one \afile\a"
"   argument, otherwise exactly one file must be specified.]"
"[+?Partitions are induced in a three step process:]{"
"   [+(1)?If \b--partition\b is not specified then a"
"   subset of columns are filtered from the training data for"
"   partitioning. This filtering usually gathers high frequency columns"
"   from the data. A high frequency column has data with a high rate of"
"   change between rows. High frequency cutoff rates, specified by"
"   \b--high\b, typically range"
"   from 10% to 50%, depending on the data. \bpzip\b(1) compresses"
"   low frequency columns much more efficiently than \bgzip\b(1).]"
"   [+(2)?A heuristic search determines an initial ordering of high"
"   frequency columns to present to step (3). An optimal solution for"
"   both ordering and partitioning is NP-complete.]"
"   [+(3)?The optimal partition for the ordering from step (2) is"
"   determined by a dynamic program that computes the compressed size"
"   for all partitions that preserve order."
"   Alternative greedy methods are also available that work more"
"   quickly but do not guarantee optimality.]"
"}"
"[+?See \bpzip\b(1) for a detailed description of file partitions and"
"   column frequencies. \bpin\b can run for a long time on some data"
"   (e.g., 10 minutes on a 400Mhz Pentium with \b--high 80 --window 4M\b)."
"   Use \b--verbose\b possibly with \b--test=010\b to monitor progress.]"

"[b:bzip?Use \bbzip\b(1) compression instead of the default \bgzip\b(1)."
"   \abzip\a is not fully supported by \apzip\a, pending further"
"   investigation.]"
"[c:cache?Generate some information on \afile\a that can be reused"
"   on another invocation; this information is saved in \bpin-\b\abase\a,"
"   where \abase\a is the base name (no directory) of \afile\a. Saved"
"   information includes column frequencies and singleton and pairwise"
"   column \agzip\a rates.]"
"[g:group?Sets the maximum number of columns in any partition group. Lower"
"   values speed up the the dynamic program but also may produce"
"   sub-optimal solutions.]#[columns:=row-size]"
"[h:high?Select this number of columns with the highest frequencies for the"
"   columns in the partition. If \acolumns\a is followed by `%' then"
"   columns with frequencies larger than this percentage are"
"   selected.]:[columns:=10%]"
"[l:level?Sets the \agzip\a compression level to \alevel\a. Levels range"
"   from 1 (fastest, worst compression) to 9 (slowest, best compression).]#"
"       [level:=6]"
"[m:maxhigh?Exit with exit code 3 if the number of high frequency columns"
"   exceeds \bmaxhigh\b. If \bmaxhigh\b is followed by `%' then the limit"
"   is \bmaxhigh\b percent of the total number of columns.]:[maxhigh:=40%]"
"[o:sort?Sort the window data by row before inducing the partition.]"
"[p:partition?Specifies the data row size and the high frequency column"
"   partition groups and permutation. The partition file is a sequence"
"   of lines. Comments start with # and continue to the end of the line."
"   The first non-comment line specifies the optional name string"
"   in \"...\". The next non-comment line specifies the row size."
"   The remaining lines operate on column offset ranges of the form:"
"   \abegin\a[-\aend\a]] where \abegin\a is the beginning column offset"
"   (starting at 0), and \aend\a is the ending column offset for an"
"   inclusive range. The operators are:]:[file]{"
"       [+range [...]]?places all columns in the specified \arange\a"
"           list in the same high frequency partition group."
"           Each high frequency partition group is processed as"
"           a separate block by the underlying compressor"
"           (\bgzip\b(1) by default).]"
"       [+range='value'?specifies that each column in \arange\a"
"           has the fixed character value \avalue\a. C-style"
"           character escapes are valid for \avalue\a.]"
"}"
"[r:row?Specifies the input row size (number of byte columns). The row size"
"   is determined by sampling the input if not specified.]#[row-size]"
"[v:verbose?List partition search details on the standard error.]"
"[w:window?Limit the number of training data rows used to induce the"
"   partition. The window size may be decreased to accomodate an"
"   integral number of complete rows.]#[window-size:=4M]"
"[O:optimize?Choose the optimization (partitioning) method for step (3)"
"   above. The methods"
"   are:]:[method:=\foptimize_default\f]{\foptimize_methods\f}"
"[Q:regress?Generate output for regression testing, such that identical"
"   invocations with identical input files will generate the same output.]"
"[R:reorder?Choose the reordering method for step (2) above. The methods"
"   are:]:[method:=\freorder_default\f]{\freorder_methods\f}"
"[S:size?Ignore \b--row\b, determine the fixed record size based on a window"
"   of sampled data, print it on the standard output, and exit. If more"
"   than one \afile\a is specified then the record size and name are"
"   printed for each file. If the sample is insufficient, or if"
"   \b--verify\b is specified, then all of the data read to determine"
"   the row size. A \b0\b size means the record size could not be"
"   determined.]"
"[T:test?Enable implementation-specific tests and tracing.]#[test-mask]{"
"   [+0x0010?Enable reorder keep trace.]"
"   [+0x0020?Enable reorder skip/cost trace.]"
"   [+0x0040?Enable reorder permutation trace.]"
"   [+0x0080?Enable reorder level 2 merge prune.]"
"   [+0x0100?Disable reorder merge prune.]"
"   [+0x0200?Partition using initial tsp cycles.]"
"}"
"[V:verify?Verify \b--size\b by reading all data instead of the window sample.]"
"[X:prefix?Uncompressed data contains a prefix that is defined by \acount\a"
"   and an optional \aterminator\a. This data is not \bpzip\b compressed."
"   \aterminator\a may be one of:]:[count[*terminator]]]{"
"       [+\aomitted\a?\acount\a bytes.]"
"       [+L?\acount\a \bnewline\b terminated records.]"
"       [+'\achar\a'?\acount\a \achar\a terminated records.]"
"}"

"\n"
"\nfile ...\n"
"\n"
"[+SEE ALSO?\bbzip\b(1), \bgzip\b(1), \bpop\b(1), \bpzip\b(1), \bpzip\b(3)]"
;

#include <ast.h>
#include <ctype.h>
#include <error.h>
#include <ls.h>
#include <pzip.h>
#include <zlib.h>
#include <bzlib.h>

#include "FEATURE/local"

#if _hdr_tsp && _lib_tspopen
#include <tsp.h>
#endif

#define OP_optimize 0x01
#define OP_reorder  0x02
#define OP_size     0x04
#define OP_verify   0x10

#define minof(x,y)  ((x)<(y)?(x):(y))

typedef struct
{
    size_t      elements;
    size_t      size;
    size_t      rate;
    size_t      hit;
    size_t      skip;
    size_t*     member;
} Part_t;

typedef struct
{
    unsigned long   frequency;
    int     prev;
    int     column;
} Stats_t;

struct Optimize_method_s;

typedef void (*Optimize_f)(struct Optimize_method_s*, unsigned char*, unsigned char*, size_t**, size_t*, size_t*, size_t, size_t, size_t);

typedef struct Optimize_method_s
{
    const char* name;
    const char* description;
    Optimize_f  fun;
} Optimize_method_t;

struct Reorder_method_s;

typedef void (*Reorder_f)(struct Reorder_method_s*, unsigned char*, unsigned char*, int*, size_t, size_t);

typedef struct Reorder_method_s
{
    const char* name;
    const char* description;
    Reorder_f   fun;
} Reorder_method_t;

static struct
{
    Stats_t*    stats;
    size_t*     pam;
    size_t*     map;
    char*       input;
    char*       cachefile;
    int     bzip;
    int     cache;
    int     level;
    int     pairs;
    int     sort;
    int     test;
    int     window;
    int     verbose;
} state;

/*
 * allocate an i-length vector of size_t
 */

static size_t*
vector(int i)
{
    register size_t*    p;

    if (!(p = newof(0, size_t, i, 0)))
        error(ERROR_SYSTEM|3, "out of space [%d vector]", i);
    return p;
}

/*
 * allocate an i X i matrix of size_t
 */

static size_t**
matrix(int i)
{
    register size_t**   p;
    register size_t*    v;
    register int        k;
    size_t          n;

    n = i * i;
    if (!(p = newof(0, size_t*, i, n * sizeof(size_t))))
        error(ERROR_SYSTEM|3, "out of space [%d X %d matrix]", i, i);
    v = (size_t*)(p + i);
    for (k = 0; k < i; k++)
    {
        p[k] = v;
        v += i;
    }
    return p;
}

/*
 * allocate an i X i partition workspace
 */

static Part_t*
partition(int i)
{
    register Part_t*    p;
    register size_t*    v;
    register int        k;
    size_t          n;

    n = i * i;
    if (!(p = newof(0, Part_t, n, n * i * sizeof(size_t))))
        error(ERROR_SYSTEM|3, "out of space [%d partition]", i);
    v = (size_t*)(p + n);
    for (k = 0; k < n; k++)
    {
        p[k].member = v;
        v += i;
    }
    return p;
}

/*
 * order frequencies hi to lo
 */

static int
byfrequency(const void* va, const void* vb)
{
    register Stats_t*   a = (Stats_t*)va;
    register Stats_t*   b = (Stats_t*)vb;

    if (a->frequency < b->frequency)
        return 1;
    if (a->frequency > b->frequency)
        return -1;
    if (a->column < b->column)
        return -1;
    if (a->column > b->column)
        return 1;
    return 0;
}

/*
 * order column lo to hi
 */

static int
bycolumn(const void* va, const void* vb)
{
    register Stats_t*   a = (Stats_t*)va;
    register Stats_t*   b = (Stats_t*)vb;

    if (a->column < b->column)
        return -1;
    if (a->column > b->column)
        return 1;
    return 0;
}

/*
 * order partition rate lo to hi
 */

static int
byrate(const void* va, const void* vb)
{
    register Part_t*    a = (Part_t*)va;
    register Part_t*    b = (Part_t*)vb;

    if (a->rate < b->rate)
        return -1;
    if (a->rate > b->rate)
        return 1;
    if (a->member < b->member)
        return -1;
    if (a->member > b->member)
        return 1;
    return 0;
}

/*
 * order row lo to hi
 */

static int
byrow(const void* va, const void* vb)
{
    return memcmp(va, vb, state.sort);
}

/*
 * dump one partition line with label to sp
 */

static void
dumppart(Sfio_t* sp, const char* label, register Part_t* pp)
{
    register int    i;

    sfprintf(sp, "reorder %s %5u %5u :", label, pp->rate, pp->size);
    for (i = 0; i < pp->elements; i++)
        sfprintf(sp, " %3u", state.map[pp->member[i]]);
    sfprintf(sp, "\n");
}

/*
 * return the compressed size of buffer b of size bytes
 */

static size_t
zip(unsigned char* b, size_t size)
{
    unsigned long       used;
    unsigned int        dest;
    size_t          safe;

    static unsigned char*   buf;
    static int      bufsize;

    safe = 4 * roundof(size, 8 * 1024);
    if (bufsize < safe)
    {
        bufsize = safe;
        if (!(buf = newof(buf, unsigned char, bufsize, 0)))
            error(ERROR_SYSTEM|3, "out of space [buf]");
    }
    if (state.bzip)
    {
        dest = bufsize;
        if (bzBuffToBuffCompress((char*)buf, &dest, (char*)b, size, state.level, 0, 30) != BZ_OK)
            error(3, "bzip compress failed");
        used = dest;
    }
    else
    {
        used = bufsize;
        if (compress2(buf, &used, b, size, state.level) != Z_OK)
            error(3, "gzip compress failed");
    }
    return used;
}

/*
 * copy col i through col j from s into t
 * and return the compressed size
 */

static size_t
field(unsigned char* t, register unsigned char* s, int i, int j, register int row, size_t tot)
{
    register unsigned char* b;
    register unsigned char* e;
    register int        n;

    b = t;
    e = s + tot;
    n = j - i + 1;
    for (s += i; s < e; s += row)
    {
        memcpy(b, s, n);
        b += n;
    }
    return zip(t, b - t);
}

/*
 * copy col i and col j from s into t
 * and return the compressed size
 */

static size_t
pair(unsigned char* t, register unsigned char* s, int i, register int j, size_t row, size_t tot)
{
    register unsigned char* b;
    register unsigned char* e;

    b = t;
    e = s + tot;
    j -= i;
    for (s += i; s < e; s += row)
    {
        *b++ = *s;
            *b++ = *(s + j);
    }
    return zip(t, b - t);
}

/*
 * return the compressed size for partition pp
 */

static size_t
part(unsigned char* t, register unsigned char* s, register Part_t* pp, size_t row, size_t tot)
{
    register unsigned char* b;
    register unsigned char* e;
    register int        i;

    b = t;
    e = s + tot;
    for (; s < e; s += row)
        for (i = 0; i < pp->elements; i++)
            *b++ = s[pp->member[i]];
    return zip(t, b - t);
}

/*
 * merge col i with part pp and keep the best compression in np
 * (np+1) is used as tmp workspace
 * return 1 if the best merge is better than i and pp separately
 */

static int
merge(unsigned char* t, unsigned char* s, int i, Part_t* pp, register Part_t* np, size_t** siz, size_t row, size_t tot)
{
    register int        j;
    register int        k;
    register Part_t*    bp;
    size_t          x;
    size_t          z;

    k = 0;
    np->member[k++] = i;
    for (j = 0; j < pp->elements; j++)
        np->member[k++] = pp->member[j];
    np->elements = k;
    bp = 0;
    x = siz[i][i] + pp->size;
    for (i = 0;;)
    {
        z = part(t, s, np, row, tot);
        if (z < x)
        {
            bp = np + 1;
            bp->size = x = z;
            for (j = 0; j < k; j++)
                bp->member[j] = np->member[j];
        }
        if (++i >= k)
            break;
        j = np->member[i];
        np->member[i] = np->member[i - 1];
        np->member[i - 1] = j;
    }
    if (bp)
    {
        np->size = bp->size;
        np->rate = np->size / k;
        for (j = 0; j < k; j++)
            np->member[j] = bp->member[j];
        if (state.test & 0x0040)
            dumppart(sfstderr, "merg", np);
        return 1;
    }
    return 0;
}

/*
 * filter out the high frequency columns from dat
 */

static size_t
filter(Sfio_t* ip, unsigned char** bufp, unsigned char** datp, Pz_t* pz, int high, int maxhigh, size_t row, size_t tot)
{
    register int    i;
    register int    j;
    char*       q;
    unsigned char*  s;
    unsigned char*  t;
    unsigned char*  dat;
    unsigned char*  buf;
    int     c;
    size_t      rows;
    size_t      n;
    ssize_t     r;
    unsigned long   freq;
    Pzpart_t*   pp;
    Sfio_t*     sp;

    if (pz || high)
    {
        buf = *bufp;
        dat = *datp;
        rows = tot / row;
        if (pz)
        {
            pp = pz->part;
            high = pp->nmap;
            if (!(state.map = newof(state.map, size_t, high, 0)))
                error(ERROR_SYSTEM|3, "out of space [map]");
            memcpy(state.map, pp->map, high * sizeof(state.map[0]));
            pzclose(pz);
        }
        else
        {
            if (!(state.map = oldof(state.map, size_t, row, 0)))
                error(ERROR_SYSTEM|3, "out of space [map]");
            if (!(state.stats = oldof(state.stats, Stats_t, row, 0)))
                error(ERROR_SYSTEM|3, "out of space [stats]");
            memset(state.stats, 0, row * sizeof(state.stats[0]));
            if (state.cache && (sp = sfopen(NiL, state.cachefile, "r")))
            {
                int     column;
                unsigned long   frequency;

                if (state.verbose)
                    sfprintf(sfstderr, "filter using %s for frequencies\n", state.cachefile);
                error_info.file = state.cachefile;
                error_info.line = 0;
                while (q = sfgetr(sp, '\n', 1))
                {
                    error_info.line++;
                    if (*q == 'f')
                    {
                        if (sfsscanf(q + 1, "%d %d %lu", &c, &column, &frequency) != 3 || c < 0 || c >= row || column < 0 || column >= row)
                            error(3, "%s: `%s' is invalid", buf, q);
                        state.stats[c].prev = 1;
                        state.stats[c].column = column;
                        state.stats[c].frequency = frequency;
                    }
                }
                for (j = 0; j < row; j++)
                    if (!state.stats[j].prev)
                        error(3, "invalid cache -- column %d frequency omitted", j);
                sfclose(sp);
                error_info.file = 0;
                error_info.line = 0;
            }
            else
            {
                if (state.verbose)
                {
                    if (high < 0)
                        sfprintf(sfstderr, "filter top %d%% high frequency columns\n", -high);
                    else if (high > 0)
                        sfprintf(sfstderr, "filter top %d high frequency columns\n", high);
                }
                s = dat;
                for (i = 0; i < rows; i++)
                    for (j = 0; j < row; j++)
                        if (state.stats[j].prev != (c = *s++))
                        {
                            state.stats[j].prev = c;
                            state.stats[j].frequency++;
                        }
                n = rows;
                while ((r = sfread(ip, s = buf, state.window)) > 0)
                {
                    r /= row;
                    if (state.sort)
                        qsort(s, r, row, byrow);
                    for (i = 0; i < r; i++)
                        for (j = 0; j < row; j++)
                            if (state.stats[j].prev != (c = *s++))
                            {
                                state.stats[j].prev = c;
                                state.stats[j].frequency++;
                            }
                    if ((n += r) > 10000000L)
                    {
                        /*
                         * that'll do pig
                         */

                        r = 0;
                        break;
                    }
                }
                if (r < 0)
                    error(ERROR_SYSTEM|3, "data read error");
                else if (r)
                    error(1, "last record incomplete");
                for (j = 0; j < row; j++)
                    state.stats[j].column = j;
                qsort(state.stats, row, sizeof(Stats_t), byfrequency);
                if (state.cache)
                {
                    error_info.file = state.cachefile;
                    error_info.line = 0;
                    if (!(sp = sfopen(NiL, state.cachefile, "w")))
                        error(ERROR_SYSTEM|1, "cannot write cache");
                    else
                    {
                        sfprintf(sp, "# %s cache for %s\n", error_info.id, state.input);
                        for (j = 0; j < row; j++)
                        {
                            sfprintf(sp, "f %d %d %lu\n", j, state.stats[j].column, state.stats[j].frequency);
                            error_info.line++;
                        }
                        if (sfclose(sp))
                            error(ERROR_SYSTEM|3, "cache write error");
                    }
                    error_info.file = 0;
                    error_info.line = 0;
                }
                if (state.verbose)
                    sfprintf(sfstderr, "filter done -- %u rows\n", n);
            }
            if (high < 0)
            {
                high = -high;
                freq = state.stats[0].frequency * high / 100;
                for (j = 0; j < row; j++)
                    if (state.stats[j].frequency <= freq)
                        break;
                high = j;
                if (state.verbose)
                    sfprintf(sfstderr, "%d high frequency column%s out of %d\n", high, high == 1 ? "" : "s", row);
            }
            if (maxhigh < 0)
            {
                maxhigh = -maxhigh;
                if (high > (row * maxhigh / 100))
                    error(6, "high frequency count %d exceeds %d%% of %d", high, maxhigh, row);
            }
            else if (maxhigh > 0)
            {
                if (high > maxhigh)
                    error(6, "high frequency count %d exceeds %d", high, maxhigh);
            }

            /*
             * cut the original data layout some slack
             * by using the original column order
             */

            qsort(state.stats, high, sizeof(Stats_t), bycolumn);
            for (j = 0; j < high; j++)
                state.map[j] = state.stats[j].column;
            state.pairs = row;
            if (!(state.pam = newof(state.pam, size_t, row, 0)))
                error(ERROR_SYSTEM|3, "out of space [pam]");
            for (j = 0; j < row; j++)
                state.pam[j] = row;
            for (j = 0; j < high; j++)
                state.pam[state.map[j]] = j;
        }
        s = dat;
        t = buf;
        for (i = 0; i < rows; i++)
        {
            for (j = 0; j < high; j++)
                *t++ = s[state.map[j]];
            s += row;
        }
        row = high;
        *bufp = dat;
        *datp = buf;
    }
    else
    {
        if (!(state.map = newof(0, size_t, row, 0)))
            error(ERROR_SYSTEM|3, "out of space [map]");
        for (i = 0; i < row; i++)
            state.map[i] = i;
    }
    return row;
}

/*
 * permute state.map and dat according to a new order
 */

static void
permute(unsigned char* buf, unsigned char* dat, size_t* ord, size_t row, size_t tot)
{
    register int    i;
    size_t*     tmap;
    unsigned char*  end;

    tmap = vector(row);
    for (i = 0; i < row; i++)
        tmap[i] = state.map[i];
    for (i = 0; i < row; i++)
        state.map[i] = tmap[ord[i]];
    for (end = dat + tot; dat < end; dat += row)
    {
        memcpy(buf, dat, row);
        for (i = 0; i < row; i++)
            dat[i] = buf[ord[i]];
    }
    free(tmap);
}

/*
 * stuff the partion group labels in lab
 */

static int
solution(int* lab, int label, size_t* sol, int beg, int end)
{
    int i;

    i = sol[end];
    if (i >= beg)
        label = solution(lab, label, sol, beg, i);
    for (label++, beg = i + 1; beg <= end; beg++)
        lab[beg] = label;
    return label;
}

/*
 * list the final partition on sp
 */

static void
list(Sfio_t* sp, int* lab, size_t row)
{
    register int    i;
    register int    j;
    register int    g;

    g = -1;
    for (i = 0; i < row; i++)
    {
        if (g != lab[i])
        {
            g = lab[i];
            sfprintf(sp, "\n");
        }
        else
            sfprintf(sp, " ");
        for (j = i + 1; j < row && lab[j] == g && state.map[j] == (state.map[j - 1] + 1); j++);
        sfprintf(sp, "%d", state.map[i]);
        if (j > (i + 2))
        {
            i = j - 1;
            sfprintf(sp, "-%d", state.map[i]);
        }
    }
    sfprintf(sp, "\n");
}

/*
 * search for a good initial ordering
 */

static void
reorder_heuristic(Reorder_method_t* method, unsigned char* buf, unsigned char* dat, int* lab, size_t row, size_t tot)
{
    register int        i;
    register int        j;
    register int        k;
    register Part_t*    cp;
    register Part_t*    np;
    register Part_t*    tp;
    register Part_t*    bp;
    Part_t*         fp;
    size_t*         hit;
    size_t*         ord;
    ssize_t*        cst;
    char*           s;
    int         ii;
    int         jj;
    size_t          y;
    size_t          z;
    ssize_t         cy;
    Part_t*         cur;
    Part_t*         nxt;
    Part_t*         fin;
    size_t**        siz;
    Sfio_t*         sp;

    siz = matrix(row);
    cur = partition(row);
    nxt = partition(row);
    fin = partition(row);
    hit = vector(row);
    ord = vector(row);
    cst = (ssize_t*)vector(row);

    /*
     * fill in the pair compression size matrix
     * and the initial partition candidate list
     */

    k = 2;
    cp = cur;

    /*
     * check the cache
     */

    if (state.cache && (sp = sfopen(NiL, state.cachefile, "r")))
    {
        if (state.verbose)
            sfprintf(sfstderr, "reorder using %s for pair compress sizes\n", state.cachefile);
        error_info.file = state.cachefile;
        error_info.line = 0;
        while (s = sfgetr(sp, '\n', 1))
        {
            error_info.line++;
            if (*s == 'p')
            {
                if (sfsscanf(s + 1, "%d %d %I*u", &ii, &jj, sizeof(z), &z) != 3 || ii < 0 || jj < 0 || state.pairs && (ii >= state.pairs || jj >= state.pairs))
                    error(3, "%s: `%s' is invalid", buf, s);
                if (state.pairs)
                {
                    ii = state.pam[ii];
                    jj = state.pam[jj];
                }
                if (ii >= row || jj >= row)
                    continue;
                siz[ii][jj] = z;
                if (ii == jj)
                    hit[ii] = 1;
                else
                {
                    hit[ii] = hit[jj] = k;
                    cp->size = siz[ii][jj];
                    cp->rate = cp->size / 2;
                    cp->elements = 2;
                    cp->member[0] = ii;
                    cp->member[1] = jj;
                    cp++;
                }
            }
        }
        sfclose(sp);
        error_info.file = 0;
        error_info.line = 0;
    }
    sp = 0;
    for (i = 0; i < row; i++)
        if (!hit[i])
        {
            hit[i] = 1;
            if (state.cache && !error_info.file)
            {
                error_info.file = state.cachefile;
                error_info.line = 0;
                if (!(sp = sfopen(NiL, error_info.file, "a")))
                    error(ERROR_SYSTEM|1, "cannot update cache");
                else if (!sfseek(sp, (Sfoff_t)0, SEEK_END))
                {
                    sfprintf(sp, "# %s cache for %s\n", error_info.id, state.input);
                    error_info.line++;
                }
            }
            siz[i][i] = field(buf, dat, i, i, row, tot);
            if (sp)
            {
                sfprintf(sp, "p %d %d %I*u\n", state.map[i], state.map[i], sizeof(siz[i][i]), siz[i][i]);
                error_info.line++;
            }
            if (state.verbose)
                sfprintf(sfstderr, "reorder pairs for %d [%d]\n", i, state.map[i]);
            for (j = 0; j < i; j++)
            {
                z = pair(buf, dat, i, j, row, tot);
                y = pair(buf, dat, j, i, row, tot);
                if (z <= y)
                {
                    ii = i;
                    jj = j;
                }
                else
                {
                    z = y;
                    ii = j;
                    jj = i;
                }
                if (z < (siz[i][i] + siz[j][j]))
                {
                    hit[i] = hit[j] = k;
                    cp->size = siz[ii][jj] = z;
                    cp->rate = cp->size / 2;
                    cp->elements = 2;
                    cp->member[0] = ii;
                    cp->member[1] = jj;
                    cp++;
                    if (sp)
                    {
                        sfprintf(sp, "p %d %d %I*u\n", state.map[ii], state.map[jj], sizeof(z), z);
                        error_info.line++;
                    }
                }
            }
        }
    if (sp && sfclose(sp))
        error(ERROR_SYSTEM|1, "cache write error");
    error_info.file = 0;
    error_info.line = 0;

    /*
     * coalesce the partitions
     */

    fp = fin;
    for (;;)
    {
        if (state.verbose)
            sfprintf(sfstderr, "reorder part %d\n", k + 1);
        qsort(cur, cp - cur, sizeof(Part_t), byrate);
        np = nxt;
        if (k == 2)
        {
            for (i = 0; i < row; i++)
                cst[i] = state.window;
            for (tp = cur; tp < cp; tp++)
            {
                i = tp->member[0];
                j = tp->member[1];
                if (!(z = siz[i][j]))
                    z = siz[j][i];
                if (z)
                {
                    cy = z - siz[i][i];
                    if (cy < cst[j])
                        cst[j] = cy;
                    cy = z - siz[j][j];
                    if (cy < cst[i])
                        cst[i] = cy;
                }
            }
        }
        else
            for (tp = cur; tp < cp; tp++)
                for (j = 0; j < tp->elements; j++)
                {
                    ii = 0;
                    for (i = 0; i < tp->elements; i++)
                        if (i != j)
                            np->member[ii++] = tp->member[i];
                    np->elements = ii;
                    cst[tp->member[j]] = tp->size - part(buf, dat, np, row, tot);
                    if (state.test & 0x0020)
                        sfprintf(sfstderr, "reorder cost %6u %u\n", state.map[tp->member[j]], cst[tp->member[j]]);
                }
        for (tp = cur; tp < cp; tp++)
        {
            bp = 0;
            for (i = 0; i < row; i++)
            {
                if (hit[i] > k)
                    continue;
                z = 0;
                jj = 0;
                for (j = 0; j < tp->elements; j++)
                    if (hit[tp->member[j]] > k)
                        z = 2;
                    else
                    {
                        tp->member[jj++] = tp->member[j];
                        if (i == tp->member[j])
                            z = 2;
                        else if (!z && (siz[i][tp->member[j]] || siz[tp->member[j]][i]))
                            z = 1;
                    }
                tp->elements = jj;
                if (z != 1)
                    continue;
                if (!merge(buf, dat, i, tp, np, siz, row, tot))
                {
                    if (!(state.test & 0x0100) && (k > 2 || (state.test & 0x0080)))
                        for (j = 0; j < tp->elements; j++)
                            siz[i][tp->member[j]] = siz[tp->member[j]][i] = 0;
                }
                else if ((ssize_t)(np->size - tp->size) < cst[i] && (!bp || np->size < bp->size))
                {
                    bp = np + 2;
                    bp->size = np->size;
                    bp->elements = np->elements;
                    for (j = 0; j < np->elements; j++)
                        bp->member[j] = np->member[j];
                }
            }
            if (bp)
            {
                np->size = bp->size;
                np->elements = bp->elements;
                for (j = 0; j < bp->elements; j++)
                    hit[np->member[j] = bp->member[j]] = k + 1;
                if (state.test & 0x0010)
                    dumppart(sfstderr, "keep", np);
                np++;
            }
            else
                tp->skip = 1;
        }
        for (tp = cur; tp < cp; tp++)
            if (tp->skip)
            {
                if (state.test & 0x0020)
                    dumppart(sfstderr, "skip", tp);
                tp->skip = 0;
                if (k > 2 && tp->elements > 1)
                {
                    fp->hit = k;
                    fp->size = tp->size;
                    fp->elements = tp->elements;
                    for (j = 0; j < tp->elements; j++)
                        fp->member[j] = tp->member[j];
                    fp++;
                }
            }
        if (np == nxt)
            break;
        tp = nxt;
        nxt = cur;
        cur = tp;
        cp = np;
        k++;
    }

    /*
     * collect the order in ord and the partition labels in lab
     */

    if (state.verbose)
        sfprintf(sfstderr, "reorder done\n");
    j = 0;
    jj = 0;
    k = 0;
    while (fp-- > fin)
        for (i = 0; i < fp->elements; i++)
            if (hit[fp->member[i]] == fp->hit)
            {
                if (!jj)
                {
                    jj = 1;
                    j++;
                }
                hit[fp->member[i]] = 0;
                ord[k++] = fp->member[i];
                lab[fp->member[i]] = j;
            }
    for (i = 0; i < row; i++)
        if (hit[i])
        {
            ord[k++] = i;
            lab[i] = ++j;
        }

    /*
     * permute state.map and dat according to the new order
     */

    permute(buf, dat, ord, row, tot);

    /*
     * clean up
     */

    free(siz);
    free(cur);
    free(nxt);
    free(fin);
    free(hit);
    free(cst);
    free(ord);
}

/*
 * tsp ordering
 */

static void
reorder_tsp(Reorder_method_t* method, unsigned char* buf, unsigned char* dat, int* lab, size_t row, size_t tot)
{
#if TSP_VERSION
    size_t      i;
    size_t      j;
    size_t      end;
    size_t      breakat;
    size_t*     self;
    size_t**    apart;
    size_t**    together;
    Tsp_t*      tsp;
    Tsp_index_t*    tour;
    Tsp_index_t*    cycle;
    Tsp_cost_t  breakval;
    Tsp_cost_t**    cost;
    Tsp_cost_t* v;
    Tsp_disc_t  disc;

    if (!(cost = newof(0, Tsp_cost_t*, row, row * row * sizeof(Tsp_cost_t))))
        error(ERROR_SYSTEM|3, "out of space [%d X %d cost matrix]", row, row);
    v = (Tsp_cost_t*)(cost + row);
    for (i = 0; i < row; i++)
    {
        cost[i] = v;
        v += row;
    }
    self = vector(row);
    apart = matrix(row);
    together = matrix(row);

    /*
     * compute the Tsp_cost_t matrix
     */

    if (state.verbose)
        sfprintf(sfstderr, "compute the tsp cost matrix\n");
    for (i = 0; i < row; i++)
        self[i] = field(buf, dat, i, i, row, tot);
    for (i = 0; i < row; i++)
        for (j = 0; j < row; j++)
            if (i != j)
            {
                together[i][j] = pair(buf, dat, i, j, row, tot);
                apart[i][j] = self[i] + self[j];
                cost[i][j] = minof(together[i][j], apart[i][j]);
            }

    /*
     * generate a tour
     */

    if (state.verbose)
        sfprintf(sfstderr, "generate a tour\n");
    memset(&disc, 0, sizeof(disc));
    disc.version = TSP_VERSION;
    disc.errorf = errorf;
    if (!(tsp = tspopen(&disc, cost, row, TSP_DFS|(state.verbose ? TSP_VERBOSE : 0))))
        error(3, "tspopen error");
    if (!(tour = tsptour(tsp)))
        error(3, "tsptour error");

    /*
     * break tour at most expensive link; put order into self
     */

    breakat = end = row-1;
    breakval = cost[tour[end]][tour[0]];
    for (i = 0; i < end; i++)
        if (cost[tour[i]][tour[i+1]] > breakval)
        {
            breakat = i;
            breakval = cost[tour[i]][tour[i+1]];
        }
    j = 0;
    for (i = breakat + 1; i < row; i++)
        self[j++] = tour[i];
    for (i = 0; i <= breakat; i++)
        self[j++] = tour[i];

    /*
     * permute state.map and dat according to the new order
     */

    permute(buf, dat, self, row, tot);

    /*
     * partition
     */

    if (state.test & 0x0200)
    {
        /*
         * partition according to the initial tsp cycles
         */

        if (!(cycle = tspcycle(tsp)))
            error(3, "tspcycle error");
        for (i = 0; i < row; i++)
            lab[i] = cycle[self[i]];
    }
    else
    {
        /*
         * partition according to dependence along tour
         */

        for (i = j = 0; i < end; i++)
        {
            lab[i] = j;
            if (together[self[i]][self[i+1]] > apart[self[i]][self[i+1]])
                j++;
        }
        lab[i] = j;
    }

    /*
     * clean up
     */

    free(self);
    free(apart);
    free(together);
    free(cost);
    tspclose(tsp);
#else
    error(3, "%s ordering requires -l%s", method->name, method->name);
#endif
}

/*
 * dynamic program to find optimal partition based on zip sizes in val
 */

static void
optimize_dynamic(Optimize_method_t* method, unsigned char* buf, unsigned char* dat, size_t** val, size_t* cst, size_t* sol, size_t row, size_t tot, size_t grp)
{
    int i;
    int j;
    int k;
    size_t  new;

    /*
     * fill in the field compress value matrix
     */

    for (i = 0; i < row; i++)
    {
        if (grp <= 0)
            k = row;
        else if ((k = i + grp) > row)
            k = row;
        if (state.verbose)
            sfprintf(sfstderr, "%s %d..%d\n", method->name, i, k - 1);
        for (j = i; j < k; j++)
            val[i][j] = field(buf, dat, i, j, row, tot);
        for (; j < row; j++)
            val[i][j] = ~0;
    }

    /*
     * now run the dynamic program
     */

    for (i = 0; i < row; i++)
    {
        cst[i] = val[0][i];
        sol[i] = -1;
        for (j = 0; j < i; j++)
        {
            new = cst[j] + val[j+1][i];
            if (new < cst[i])
            {
                cst[i] = new;
                sol[i] = j;
            }
        }
    }
}

/*
 * greedy algorithm to approximate optimal partition based on zip sizes in val
 */

static void
optimize_greedy(Optimize_method_t* method, unsigned char* buf, unsigned char* dat, size_t** val, size_t* cst, size_t* sol, size_t row, size_t tot, size_t grp)
{
    int i;
    size_t  expand;
    size_t  new;

    cst[0] = val[0][0] = field(buf, dat, 0, 0, row, tot);
    sol[0] = -1;
    for (i = 1; i < row; i++)
    {
        val[sol[i-1] + 1][i] = field(buf, dat, sol[i-1]+1, i, row, tot);
        val[i][i] = field(buf, dat, i, i, row, tot);
        expand = val[sol[i-1] + 1][i];
        new = cst[i-1] + val[i][i];
        if (new < expand) {
            cst[i] = val[i][i];
            sol[i] = i-1;
        }
        else {
            cst[i] = expand;
            sol[i] = sol[i-1];
        }
    }
}

/*
 * transitive greedy algorithm to approximate optimal partition
 * based on zip sizes in val
 */

static void
optimize_transitive(Optimize_method_t* method, unsigned char* buf, unsigned char* dat, size_t** val, size_t* cst, size_t* sol, size_t row, size_t tot, size_t grp)
{
    int i;
    size_t  expand;

    cst[0] = val[0][0] = field(buf, dat, 0, 0, row, tot);
    sol[0] = -1;
    for (i = 1; i < row; i++)
    {
        val[i][i] = field(buf, dat, i, i, row, tot);
        expand = field(buf, dat, i-1, i, row, tot);
        if (expand <= val[i-1][i-1] + val[i][i])
        {
            cst[i] = expand;
            sol[i] = sol[i-1];
        }
        else
        {
            cst[i] = val[i][i];
            sol[i] = i-1;
        }
    }
}

/*
 * optimization (partitioning) method table
 * the first entry is the default
 */

static Optimize_method_t    optimize_methods[] =
{
    {
        "dynamic",
        "dynamic program optimal partition",
        optimize_dynamic
    },
    {
        "greedy",
        "greedy approximation partition",
        optimize_greedy
    },
    {
        "none",
        "no partition",
        0
    },
    {
        "transitive",
        "transitive greedy approximation partition",
        optimize_transitive
    },
};

/*
 * reorder method table
 * the first entry is the default
 */

static Reorder_method_t     reorder_methods[] =
{
    {
        "heuristic",
        "heuristic reorder",
        reorder_heuristic
    },
    {
        "none",
        "no reordering",
        0
    },
    {
        "tsp",
        "tsp reordering",
        reorder_tsp
    },
};

/*
 * optget() info discipline function
 */

static int
optinfo(Opt_t* op, Sfio_t* sp, const char* s, Optdisc_t* dp)
{
    register int    i;
    register int    n;

    n = 0;
    if (streq(s, "optimize_default"))
        n += sfprintf(sp, "%s", optimize_methods[0].name);
    else if (streq(s, "optimize_methods"))
        for (i = 0; i < elementsof(optimize_methods); i++)
            n += sfprintf(sp, "[+%s?%s]", optimize_methods[i].name, optimize_methods[i].description);
    else if (streq(s, "reorder_default"))
        n += sfprintf(sp, "%s", reorder_methods[0].name);
    else if (streq(s, "reorder_methods"))
        for (i = 0; i < elementsof(reorder_methods); i++)
            n += sfprintf(sp, "[+%s?%s]", reorder_methods[i].name, reorder_methods[i].description);
    return n;
}

int
main(int argc, char** argv)
{
    unsigned char*  dat;
    unsigned char*  buf;
    char*       s;
    int     i;
    ssize_t     n;
    ssize_t     m;
    size_t      rows;
    size_t      tot;
    size_t      rec;
    size_t      win;
    int*        lab;
    Sfio_t*     ip;
    Sfio_t*     dp;
    Pzdisc_t    disc;
    struct stat is;
    struct stat cs;
    Optdisc_t   optdisc;

    int     op = 0;
    int     high = -10;
    int     maxhigh = -40;
    int     maxgrp = 0;
    size_t      row = 0;
    Pz_t*       pz = 0;
    char*       partition = 0;

    Optimize_method_t*  optimize_method = &optimize_methods[0];
    Reorder_method_t*   reorder_method = &reorder_methods[0];

    error_info.id = "pin";
    optinit(&optdisc, optinfo);
    state.level = 6;
    state.window = PZ_WINDOW;
    if (!(dp = sfstropen()))
        error(ERROR_SYSTEM|3, "out of space");
    for (;;)
    {
        switch (optget(argv, usage))
        {
        case 'b':
            state.bzip = 1;
            continue;
        case 'c':
            state.cache = opt_info.num;
            continue;
        case 'g':
            maxgrp = opt_info.num;
            continue;
        case 'h':
            high = strtol(opt_info.arg, &s, 0);
            if (*s == '%')
            {
                s++;
                high = -high;
            }
            if (*s)
                error(3, "%s: %s: invalid number", opt_info.name, opt_info.arg);
            continue;
        case 'l':
            state.level = opt_info.num;
            continue;
        case 'm':
            maxhigh = strtol(opt_info.arg, &s, 0);
            if (*s == '%')
            {
                s++;
                maxhigh = -maxhigh;
            }
            if (*s)
                error(3, "%s: %s: invalid number", opt_info.name, opt_info.arg);
            continue;
        case 'o':
            state.sort = 1;
            continue;
        case 'p':
            partition = opt_info.arg;
            continue;
        case 'r':
            row = opt_info.num;
            continue;
        case 'v':
            state.verbose = 1;
            continue;
        case 'w':
            state.window = opt_info.num;
            continue;
        case 'O':
            n = strlen(opt_info.arg);
            for (i = 0;; i++)
            {
                if (i >= elementsof(optimize_methods))
                {
                    error(2, "%s: unknown partition method", opt_info.arg);
                    break;
                }
                if (strneq(opt_info.arg, optimize_methods[i].name, n))
                {
                    optimize_method = &optimize_methods[i];
                    break;
                }
            }
            continue;
        case 'Q':
            sfprintf(dp, "regress\n");
            continue;
        case 'R':
            n = strlen(opt_info.arg);
            for (i = 0;; i++)
            {
                if (i >= elementsof(reorder_methods))
                {
                    error(2, "%s: unknown reorder method", opt_info.arg);
                    break;
                }
                if (strneq(opt_info.arg, reorder_methods[i].name, n))
                {
                    reorder_method = &reorder_methods[i];
                    break;
                }
            }
            continue;
        case 'S':
            op |= OP_size;
            continue;
        case 'T':
            sfprintf(dp, "test=%s\n", opt_info.arg);
            continue;
        case 'V':
            op |= OP_verify;
            continue;
        case 'X':
            sfprintf(dp, "prefix=%s\n", opt_info.arg);
            continue;
        case '?':
            error(ERROR_USAGE|4, "%s", opt_info.arg);
            continue;
        case ':':
            if (!opt_info.option[0] && opt_info.name[0] == opt_info.name[1])
                sfputr(dp, &argv[opt_info.index - 1][2], '\n');
            else
                error(2, "%s", opt_info.arg);
            continue;
        }
        break;
    }
    argv += opt_info.index;
    if (error_info.errors || !(state.input = *argv++) || (n = *argv != 0) && !(op & OP_size))
        error(ERROR_USAGE|4, "%s", optusage(NiL));

    /*
     * set up the workspace
     */

    if (optimize_method->fun)
        op |= OP_optimize;
    if (reorder_method->fun)
        op |= OP_reorder;
    memset(&disc, 0, sizeof(disc));
    disc.version = PZ_VERSION;
    disc.errorf = errorf;
    disc.partition = partition;
    disc.window = 1;
    if (sftell(dp) && !(disc.options = strdup(sfstruse(dp))))
        error(ERROR_SYSTEM|3, "out of space");
    sfclose(dp);
    if (op & OP_size)
    {
        row = 0;
        dat = 0;
        buf = 0;
    }
    for (;;)
    {
        ip = 0;
        if (pz = pzopen(&disc, state.input, 0))
        {
            state.test = pz->test;
            ip = pz->io;
            if (partition)
            {
                row = pz->part->row;
                high = pz->part->nmap;
            }
            else
            {
                m = pz->part->row;
                if (!row)
                {
                    if (m > 0)
                    {
                        row = m;
                        if (state.verbose && !(op & OP_size))
                            sfprintf(sfstderr, "row size %I*u\n", sizeof(row), row);
                    }
                    else if (!(op & OP_size))
                        error(1, "could not determine row size");
                }
                else if (m > 0 && (row < m || row > m && (row % m)))
                    error(1, "row size %I*d may be invalid -- try %I*u next time", sizeof(row), row, sizeof(m), m);
                pz->io = 0;
                pzclose(pz);
                pz = 0;
            }
        }
        else if (!(op & OP_size))
            return 1;
        if (!(op & OP_size))
            break;
        if (row && (op & OP_verify))
        {
            win = (state.window / row) * row;
            if (!(dat = newof(dat, unsigned char, win, 0)))
                error(ERROR_SYSTEM|3, "out of space [dat]");
            if (!(buf = newof(buf, unsigned char, win, 0)))
                error(ERROR_SYSTEM|3, "out of space [buf]");
            if ((m = sfread(ip, dat, win)) < 0)
                m = 0;
            rows = m / row;
            if (filter(ip, &buf, &dat, pz, high, maxhigh, row, row * rows) == row)
                row = 0;
        }
        if (!n)
        {
            sfprintf(sfstdout, "%I*u\n", sizeof(row), row);
            return 0;
        }
        sfprintf(sfstdout, "%8I*u %s\n", sizeof(row), row, state.input);
        if (pz)
            pzclose(pz);
        else if (ip)
            sfclose(ip);
        if (!(state.input = *argv++))
            return 0;
        row = 0;
    }
    if (!(rec = row))
        error(3, "-r row-size is required");
    if (high > (int)row)
        error(3, "-h col-count must be <= -r row-size");
    state.window = (state.window / row) * row;
    if (!(dat = newof(0, unsigned char, state.window, 0)))
        error(ERROR_SYSTEM|3, "out of space [dat]");
    if (!(buf = newof(0, unsigned char, state.window, 0)))
        error(ERROR_SYSTEM|3, "out of space [buf]");
    if ((n = sfread(ip, dat, state.window)) <= 0)
        error(3, "input empty");
    rows = n / row;
    state.window = rows * row;
    if (state.sort)
    {
        state.sort = row;
        qsort(dat, rows, row, byrow);
    }

    /*
     * set up the cache file
     */

    if (stat(state.input, &is))
        error(ERROR_SYSTEM|3, "%s: cannot stat", state.input);
    if (s = strrchr(state.input, '/'))
        s++;
    else
        s = state.input;
    if (state.cache)
    {
        if (!(state.cachefile = strdup(sfprints("%s-%s", error_info.id, s))))
            error(ERROR_SYSTEM|3, "out of space");
        if (!stat(state.cachefile, &cs) && cs.st_mtime < is.st_mtime && remove(state.cachefile))
            error(ERROR_SYSTEM|3, "%s: cannot update intermediate data cache", state.cachefile);
    }

    /*
     * filter out the high frequency columns
     */

    if (row = filter(ip, &buf, &dat, pz, high, maxhigh, row, row * rows))
    {
        tot = row * rows;
        if (!(lab = newof(0, int, row, 0)))
            error(ERROR_SYSTEM|3, "out of space [lab]");
        if (pz && !(op & (OP_reorder|OP_optimize)))
            for (n = 0; n < pz->part->nmap; n++)
                lab[n] = pz->part->lab[n];
    }
    else
        op &= ~(OP_optimize|OP_reorder);

    /*
     * determine a better ordering
     */

    if (op & OP_reorder)
        (*reorder_method->fun)(reorder_method, buf, dat, lab, row, tot);

    /*
     * generate a partition based on the ordering
     */

    if (op & OP_optimize)
    {
        size_t*     cst;
        size_t*     sol;
        size_t**    val;

        /*
         * allocate the tmp workspace
         */

        val = matrix(row);
        cst = vector(row);
        sol = vector(row);

        /*
         * partition
         */

        (*optimize_method->fun)(optimize_method, buf, dat, val, cst, sol, row, tot, maxgrp);

        /*
         * gather the optimal partition group labels in lab
         */

        solution(lab, 0, sol, 0, row - 1);
        if (state.verbose)
            sfprintf(sfstderr, "%s done\n", optimize_method->name);

        /*
         * clean up
         */

        free(val);
        free(cst);
        free(sol);
    }

    /*
     * finished -- list the partition
     */

    sfprintf(sfstdout, "# pzip partition\n");
    sfprintf(sfstdout, "# %s\n", fmtident(usage));
    if (!(op & OP_reorder))
        sfprintf(sfstdout, "# group coalescing limited to adjacent columns\n");
    if (maxgrp)
        sfprintf(sfstdout, "# group size limited to %d columns\n", maxgrp);
    sfprintf(sfstdout, "# row %d window %d compression level %d\n", rec, state.window, state.level);
    if (state.sort)
        sfprintf(sfstdout, "\nsort=1\n");
    if (high)
        sfprintf(sfstdout, "\n%d\t# high frequency %d\n", rec, row);
    else
        sfprintf(sfstdout, "\n%d\n", rec);
    list(sfstdout, lab, row);
    return 0;
}