/* column.c
********************************************************************************
* column.c -- Routines for treating a PDS COLUMN object as a series.
*
* User routines:
*	Pro_ColumnSeries()	creates a series object out of a specified
*				PDS-labeled COLUMN object.
*
* Mark Showalter & Neil Heather, PDS Ring-Moon Systems Node, March 1998
*******************************************************************************/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "profile.h"
#include "fortran.h"
#include "oal.h"

/*************************
 * Data type definitions *
 *************************/

typedef struct ZPRO_COLUMN_STRUCT {
    XPRO_CLASS	class;
    PRO_OBJECT	*label;
    ODLTREE	tabletree, columntree;
    RL_INT4	table, column, item, items, row1, row2, rows;
    RL_INT4	k1, start, skip;
    RL_FLT8	missing, invalid, offset, scale;
    RL_BOOL	isdouble;
    RL_FLT8	*(ddata[2]);
    RL_FLT4	*(sdata[2]);
    RL_INT4	buffer_rows, rowsegs[2];
} ZPRO_COLUMN;

/********************
 * Static variables *
 ********************/

static XPRO_CLASS column_class = {XPRO_COLUMN_CLASS, "PDS column", NULL};

/************************************
 * Prototypes of internal functions *
 ************************************/

static RL_FLT8 ZPro_ColumnValue RL_PROTO((RL_VOID *pointer, RL_INT4 k,
                                          RL_INT4 *flag));
static void    ZPro_FreeColumn  RL_PROTO((RL_VOID *pointer));
static void    ZPro_PrintColumn RL_PROTO((RL_VOID *pointer));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_ColumnSeries (column.c)
*$ Abstract:
*	This routine creates and returns a column series object.  A column
*	series returns the values from a column of data found in a series or
*	table of a PDS-labeled data file.
*$ Keywords:
*	PROFILE, SERIES
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	PRO_OBJECT	*Pro_ColumnSeries(label, ntable, ncolumn, nitem,
*				row1, row2, k1, buffersize, usedouble)
*	PRO_OBJECT	*label;
*	RL_INT4		ntable, ncolumn, nitem, row1, row2, k1, buffersize;
*	RL_BOOL		usedouble;
*$ Inputs:
*	label		pointer to a PDS label object.
*	ntable		table index in label.
*	ncolumn		column index in label.
*	nitem		item index [1...n] or 0 to use all items.
*	row1		first row to include (numbered from 1).
*	row2		last row to include (numbered from 1; 0 for all).
*	k1		starting index for series.
*	buffersize	minimum number of samples to keep in internal buffer;
*			0 to keep entire column in memory.
*	usedouble	TRUE to save double precision values internally;
*			FALSE for single precision.
*$ Outputs:
*	none
*$ Returns:
*	pointer to a new column series object, or NULL on non-fatal error.
*$ Detailed_description:
*	This routine creates and returns a column series object.  A column
*	series returns the values from a column of data found in a series or
*	table of a PDS-labeled data file.
*
*	For PDS SERIES objects, the SAMPLING_PARAMETER quantities in the label
*	define the sampling of the series.  For PDS TABLE objects, the series
*	sampling parameter is assumed to be the row index, starting with 1.
*
*	Note: The buffersize parameter enables the user to have access to a
*	finite subset of the array at any given time without the need to retain
*	the entire column in memory.  Whenever an attempt is made to access a
*	sample not currently in memory, the needed segment of the data file is
*	read into memory.  The buffersize parameter should be no smaller than
*	the difference between the first and last series samples needed at a
*	given time.  For example, when convolving the data with a filter
*	function, the buffersize should be greater than or equal to the full
*	width of the filter used.
*$ External_references:
*	Profile toolkit, Object Access library
*$ Side_effects:
*	Memory is allocated.  A link to the corresponding label object is
*	created.
*$ Examples:
*	This snippet of code prints out the first three samples in the second
*	column of the first table in a PDS-labeled data file.
*
*	label = Pro_OpenLabel("table.tab");
*	ntable = 1;
*	ncolumn = 2;
*	nitem = 0;
*	k1 = 0;
*	column = Pro_ColumnSeries(label, ntable, ncolumn, nitem,
*				1, 0, k1, 0, TRUE);
*
*	printf("%g %g %g\n", Pro_SeriesValue(column,k1),
*	                     Pro_SeriesValue(column,k1+1),
*	                     Pro_SeriesValue(column,k1+2));
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if label is NULL or is not a label object.
*	PRO_DOMAIN_ERROR	it ntable, ncolumn or nitem is out of range.
*	PRO_MEMORY_ERROR	on memory allocation failure.
*	PRO_SETUP_FAILURE	on any OA Library error.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

PRO_OBJECT	*Pro_ColumnSeries(label, ntable, ncolumn, nitem, row1, row2, k1,
			buffersize, usedouble)
PRO_OBJECT	*label;
RL_INT4		ntable, ncolumn, nitem, row1, row2, k1, buffersize;
RL_BOOL		usedouble;
{
PRO_OBJECT	*new;
ZPRO_COLUMN	*column;
RL_INT4		k2, rows, flag, oldtype;
RL_FLT8		x1, x2, dx;
RL_CHAR		*xname, *yname;
RL_BOOL		status, olderror;
ODLTREE		nodes[1];
RL_VOID		*pointer, *tableptr, *columnptr;

    /* Allocate new column structure */
    column = (ZPRO_COLUMN *) XRL_Malloc(sizeof(ZPRO_COLUMN));
    if (column == NULL) return NULL;

    /* Initialize structure */
    column->class    = column_class;
    column->label    = label;
    column->table    = ntable;
    column->column   = ncolumn;
    column->k1       = k1;
    column->isdouble = usedouble;
    column->ddata[0] = NULL;
    column->ddata[1] = NULL;
    column->sdata[0] = NULL;
    column->sdata[1] = NULL;
    column->rowsegs[0] = 1;
    column->rowsegs[1] = 0;

    /* Get column info from label */
    status = XPro_LabelInfo(label, ntable, ncolumn, &tableptr, &columnptr,
		&(column->rows), &(column->items),
		&x1, &x2, &dx, &xname, &yname,
		&(column->missing), &(column->invalid),
		&(column->offset),  &(column->scale));
    if (!status) {
	XRL_Free((RL_VOID *) column);
	return NULL;
    }

    column->tabletree  = (ODLTREE) tableptr;
    column->columntree = (ODLTREE) columnptr;

    if (!usedouble) {
	if (column->missing != -HUGE_VAL)
		column->missing = (RL_FLT4) column->missing;
	if (column->invalid != -HUGE_VAL)
		column->invalid = (RL_FLT4) column->invalid;
    }

    /* Interpret row parameters */
    if (row1 < 1 || row1 > column->rows) {
	XPro_IDomainError("first row index", label, 1, column->rows, row1);
	XRL_Free((RL_VOID *) column);
	return NULL;
    }
    if (row2 == 0) row2 = column->rows;
    if (row2 < row1 || row2 > column->rows) {
	XPro_IDomainError("last row index", label, row1, column->rows, row2);
	XRL_Free((RL_VOID *) column);
	return NULL;
    }
    column->row1 = row1;
    column->row2 = row2;

    /* Interpret item parameter */
    if (nitem == 1 && column->items == 1) nitem = 0;

    if (nitem == 0) {
	dx /= column->items;
	k2 = k1 + (row2-row1+1) * column->items - 1;
	column->start = 0;
	column->skip  = 1;
    }
    else {
	if (nitem < 1 || nitem > column->items) {
	    XPro_IDomainError("item index", label, 1, column->items, nitem);
	    XRL_Free((RL_VOID *) column);
	    return NULL;
	}
	column->start = nitem-1;
	column->skip  = column->items;
	k2 = k1 + row2 - row1;
    }
    column->item = nitem;
    column->k1 = k1;

    /* Determine number of rows to retain in cache */
    rows = row2 - row1 + 1;
    if (buffersize <= 0)
	column->buffer_rows = rows;
    else if (column->skip > 1)
	column->buffer_rows = buffersize;
    else
	column->buffer_rows = (buffersize + column->items - 1) / column->items;

    if (column->buffer_rows > rows) column->buffer_rows = rows;

    /* Save any PRO_EVALUATION_FAILURE and record future ones */
    olderror = RL_ClearError1("PRO_EVALUATION_FAILURE");
    oldtype = RL_SetErrorType("PRO_EVALUATION_FAILURE", RL_RECORD);

    /* Fill first buffer and check for errors */
    ZPro_ColumnValue((RL_VOID *) column, column->k1, &flag);
    status = RL_ClearError1("PRO_EVALUATION_FAILURE");

    /* Raise PRO_SETUP_FAILURE if necessary */
    if (status) {
	RL_RaiseError("PRO_SETUP_FAILURE", xpro_message);
	/* Note xpro_message is a global variable and is already loaded with the
           error message by Pro_ColumnValue() */
    }

    /* Restore any PRO_EVALUATION_FAILURE */
    if (olderror) RL_RaiseError("PRO_EVALUATION_FAILURE", "");
    RL_SetErrorType("PRO_EVALUATION_FAILURE", oldtype);

    /* Free the object and return NULL on error */
    if (status) {
	ZPro_FreeColumn((RL_VOID *) column);
	return NULL;
    }

    /* Create new object */
    new = XPro_MakeSeries(k1, k2, x1, dx,
                          ZPro_ColumnValue, ZPro_FreeColumn, ZPro_PrintColumn,
                          (RL_VOID *) column);

    /* Enslave the label and check for error */
    if (new != NULL && XPro_EnslaveObject(new, label)) {
	Pro_FreeObject(new);
    }

    /* Transfer coordinate names to new object */
    Pro_RenameObject(new, 1, xname);
    Pro_RenameObject(new, 2, yname);

    return new;
}

/*
********************************************************************************
* INTERNAL FUNCTIONS
********************************************************************************
* ZPro_ColumnValue(pointer, k, flag)
*
* This is the series evaluation function for a PDS column series object.  It
* checks to see if the sample is in the buffer, and reads a new section of the
* data file if necessary.
*
* Inputs:
*	pointer		pointer to the ZPRO_COLUMN data structure.
*	k		index at which to evaluate array.
*
* Outputs:
*	*flag		0 if value returned is valid; PRO_MISSING_FLAG if it is
*			missing; PRO_INVALID_FLAG if it is invalid.
*
* Return:		value of column at given index; 0. on error.
*
* Errors:
*	PRO_EVALUATION_FAILURE	if an I/O occurs during file access.
*******************************************************************************/

static RL_FLT8	ZPro_ColumnValue(pointer, k, flag)
RL_VOID		*pointer;
RL_INT4		k, *flag;
{
ZPRO_COLUMN	*column;
RL_INT4		rowrel, rowseg, rowbase, rowmod, row1, row2, offset;
RL_FLT8		value;
RL_VOID		*array;
RL_CHAR		*verb;
OA_OBJECT	oaobject, oaconvert;

    column = (ZPRO_COLUMN *) pointer;

    /***************************************
     * Identify needed segment of the file *
     ***************************************/

    rowrel  = (k - column->k1) * column->skip/column->items;
    rowseg  = rowrel / column->buffer_rows;
    rowbase = rowseg * column->buffer_rows;
    rowmod  = rowseg & 1;

    /*
     * Explanation...
     *
     * rowrel is the desired sample's row relative to the first row used by
     *        this column object.
     * rowseg is the segment index.  This is 0 for the first set of
     *        <buffer_rows> rows in the file, 1 for the second set, etc.
     * rowbase is the first row in the segment, relative to the first row used.
     * rowmod is 0 if rowseg is even; 1 if it is odd.  This enables the routine
     *        to keep two different segments in memory simultaneously.
     */

    /*****************************************
     * Load segment into memory if necessary *
     *****************************************/

    if (rowseg != column->rowsegs[rowmod]) {
	column->rowsegs[rowmod] = rowseg;

	/* Read segment of column from PDS file */
	row1 = column->row1 + rowbase;
	row2 = row1 + column->buffer_rows - 1;
	if (row2 > column->row2) row2 = column->row2;

#ifdef DEBUG
	printf("Reading rows %d to %d\n", row1, row2);
#endif

	oaobject = OaReadSubTable(column->tabletree, (long) row1, (long) row2,
	                          &(column->columntree), 1);

	if (oaobject == NULL || (oa_errno >= 500 && oa_errno < 900)) {
	    verb = "reading";
	    goto OAFAILURE;
        }

	/* Convert the data */
	if (column->isdouble) {
	    oaconvert = OaConvertObjecttoOneType(oaobject, "double", 0, FALSE);
	} else {
            oaconvert = OaConvertObjecttoOneType(oaobject, "float",  0, FALSE);
	}

	if (oaconvert == NULL || (oa_errno >= 500 && oa_errno < 900)) {
	    verb = "converting";
	    goto OAFAILURE;
	}

	/* Export the data */
	array = OaExportObject(oaconvert);
	if (column->isdouble) {
	    XRL_Free((RL_VOID *) column->ddata[rowmod]);
	    column->ddata[rowmod] = (RL_FLT8 *) array;
	}
	else {
	    XRL_Free((RL_VOID *) column->sdata[rowmod]);
            column->sdata[rowmod] = (RL_FLT4 *) array;
	}

	if (pointer == NULL || (oa_errno >= 500 && oa_errno < 900)) {
	    verb = "exporting";
	    goto OAFAILURE;
	}

	/* Delete the unconverted object */
	(void) OaDeleteObject(oaobject);

	if (oa_errno >= 500 && oa_errno < 900) {
	    verb = "deleting";
	    goto OAFAILURE;
	}
    }

    /*******************************
     * Extract sample from segment *
     *******************************/

    offset = (k - column->k1)*column->skip - rowbase * column->items
             + column->start;
    if (column->isdouble) {
	if (column->ddata == NULL) goto INVALID; /* if I/O error occurred */
	value = (column->ddata[rowmod])[offset];
    }
    else {
	if (column->sdata == NULL) goto INVALID; /* if I/O error occurred */
	value = (RL_FLT8) (column->sdata[rowmod])[offset];
    }

    /* Check for flags */
    if (value == column->missing) goto MISSING;
    if (value == column->invalid) goto INVALID;

    /* Return valid value */
    *flag = 0;
    return value * column->scale + column->offset;

OAFAILURE:
    (void) sprintf(xpro_message, "error %s PDS column object\n\
OA Library error code %d", verb, oa_errno);
    RL_RaiseError("PRO_EVALUATION_FAILURE", xpro_message);
    goto INVALID;

INVALID:
    *flag = PRO_INVALID_FLAG;
    return 0.;

MISSING:
    *flag = PRO_MISSING_FLAG;
    return 0.;
}

/*
********************************************************************************
* ZPro_FreeColumn(pointer)
*
* This is the freeing function for a PDS column series object.
*
* Inputs:
*	pointer		pointer to the ZPRO_COLUMN data structure.
*******************************************************************************/

static void	ZPro_FreeColumn(pointer)
RL_VOID		*pointer;
{
ZPRO_COLUMN	*column;
RL_FLT8		value;

    column = (ZPRO_COLUMN *) pointer;
    XRL_Free((RL_VOID *) column->ddata[0]);
    XRL_Free((RL_VOID *) column->ddata[1]);
    XRL_Free((RL_VOID *) column->sdata[0]);
    XRL_Free((RL_VOID *) column->sdata[1]);

    XRL_Free(pointer);
}

/*
********************************************************************************
* ZPro_PrintColumn(pointer)
*
* This is the series printing function for a PDS series object.
*
* Inputs:
*	pointer		pointer to the ZPRO_COLUMN data structure.
*******************************************************************************/

static void	ZPro_PrintColumn(pointer)
RL_VOID		*pointer;
{
ZPRO_COLUMN	*column;
RL_INT4		k;

    column = (ZPRO_COLUMN *) pointer;

    /* Make sure object is not NULL */
    if (column == NULL) {
	printf("PRINT ERROR: PDS column series pointer is NULL\n");
	return;
    }

    /* Make sure object is a PDS series */
    if (column->class.id != XPRO_COLUMN_CLASS) {
	printf("PRINT ERROR: Object is not a PDS column series\n");
	return;
    }

    /* Print object info... */
    printf("\nPDS column series parameters...\n");
    printf("    label = "); XPro_PrintInfo(column->label);
    printf("    table = %d\n", column->table);
    printf("   column = %d\n", column->column);

    if (column->items == 1) {
	printf("     item = 1 of 1\n");
    } else if (column->item == 0) {
	printf("     item = all of %d\n", column->items);
    } else {
	printf("     item = %d of %d\n", column->item, column->items);
    }

    printf("     rows = %1d to %1d of %1d\n",
			column->row1, column->row2, column->rows);
    printf("precision = %s\n", column->isdouble ? "double":"single");
    if (column->missing != -HUGE_VAL)
	printf("  missing = %#g\n", column->missing);
    if (column->invalid != -HUGE_VAL)
	printf("  invalid = %#g\n", column->invalid);
    if (column->scale != 1. && column->offset != 0.)
	printf("  scaling = %#g * x + %#g\n", column->scale, column->offset);
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
*	FPro_ColumnSeries (column.c)
*$ Abstract:
*	This routine creates and returns a column series object.  A column
*	series returns the values from a column of data found in a series or
*	table of a PDS-labeled data file.
*$ Keywords:
*	PROFILE, SERIES
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	integer*4 function FPro_ColumnSeries(label, table, column, item,
*				row1, row2, k1, buffersize, usedouble)
*	integer*4	label, table, column, item, row1, row2, k1, buffersize
*	logical*4	usedouble
*$ Inputs:
*	label		FORTRAN pointer to a PDS label object.
*	ntable		table index in label.
*	ncolumn		column index in label.
*	nitem		item index [1...n] or 0 to use all items.
*	row1		first row to include (numbered from 1).
*	row2		last row to include (numbered from 1; 0 for all).
*	k1		starting index for series.
*	buffersize	minimum number of samples to keep in internal buffer;
*			0 to keep entire column in memory.
*	usedouble	TRUE to save double precision values internally;
*			FALSE for single precision.
*$ Outputs:
*	none
*$ Returns:
*	FORTRAN pointer to a new column series object, or NULL on non-fatal
*	error.
*$ Detailed_description:
*	This routine creates and returns a column series object.  A column
*	series returns the values from a column of data found in a series or
*	table of a PDS-labeled data file.
*
*	For PDS SERIES objects, the SAMPLING_PARAMETER quantities in the label
*	define the sampling of the series.  For PDS TABLE objects, the series
*	sampling parameter is assumed to be the row index, starting with 1.
*
*	Note: The buffersize parameter enables the user to have access to a
*	finite subset of the array at any given time without the need to retain
*	the entire column in memory.  Whenever an attempt is made to access a
*	sample not currently in memory, the needed segment of the data file is
*	read into memory.  The buffersize parameter should be no smaller than
*	the difference between the first and last series samples needed at a
*	given time.  For example, when convolving the data with a filter
*	function, the buffersize should be greater than or equal to the full
*	width of the filter used.
*$ External_references:
*	Profile toolkit, Object Access library
*$ Side_effects:
*	Memory is allocated.  A link to the corresponding label object is
*	created.
*$ Examples:
*	This snippet of code prints out the first three samples in the second
*	column of the first table in a PDS-labeled data file.
*
*	label = Pro_OpenLabel('table.tab')
*	ntable = 1
*	ncolumn = 2
*	nitem = 0
*	k1 = 0
*	column = FPro_ColumnSeries(label, ntable, ncolumn, nitem,
*    &				1, 0, k1, 0, .TRUE.)
*
*	write(*,*) FPro_SeriesValue(column,k1),
*    &	           FPro_SeriesValue(column,k1+1),
*    &	           FPro_SeriesValue(column,k1+2)
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if label is NULL or is not a label object.
*	PRO_DOMAIN_ERROR	it ntable, ncolumn or nitem is out of range.
*	PRO_MEMORY_ERROR	on memory allocation failure.
*	PRO_SETUP_FAILURE	on any OA Library error.
*	FORTRAN_POINTER_ERROR	if label is not a valid FORTRAN object pointer.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_INT4 FORTRAN_NAME(fpro_columnseries) (label, table, column, item, row1, row2,
				k1, buffersize, usedouble)
RL_INT4 *label, *table, *column, *item, *row1, *row2, *k1, *buffersize,
	*usedouble;
{
RL_VOID *ptr1, *ptr2;
RL_INT4 index;

    /* Look up label pointer */
    ptr1 = FORT_GetPointer(*label);
    if (ptr1 == NULL) return 0;

    /* Call function */
    ptr2 = (RL_VOID *) Pro_ColumnSeries(ptr1, *table, *column, *item, *row1,
				*row2, *k1, *buffersize, (RL_BOOL) *usedouble);
    if (ptr2 == NULL) return 0;

    index = FORT_AddPointer(ptr2);
    if (index == 0) Pro_FreeObject((PRO_OBJECT *) ptr2);

    return index;
}

/*******************************************************************************
*/