********************************************************************************
* filtered.c -- Routines to filter a statistical series objects
*
* User routines
*	Pro_FilteredStat()	creates a filtered statistical series object.
*
* A filtered statistical series object is a series where each sample is a
* weighted sum of the corresponding samples in another series.  The filter
* function is assumed to be constant, so the new series is, in effect, a
* convolution of the old series with a fixed filter.
*
* Mark Showalter, PDS Ring-Moon Systems Node, March 1998
*******************************************************************************/
#include <stdio.h>
#include "profile.h"
#include "fortran.h"

/********************
 * Type definitions *
 ********************/

typedef struct ZPRO_FILTERED_STRUCT {
    PRO_OBJECT	*stat;
    RL_INT4	df1, df2, fsize;
    RL_FLT8	*filter;
} ZPRO_FILTERED;

/********************************
 * Internal function prototypes *
 ********************************/

static void    ZPro_FilteredKRange RL_PROTO((RL_INT4 knew,
                                             RL_INT4 *k1, RL_INT4 *k2,
                                             RL_VOID *params));
static RL_FLT8 ZPro_FilteredWeight RL_PROTO((RL_INT4 knew, RL_INT4 k,
                                             RL_VOID *params));
static void    ZPro_PrintFiltered  RL_PROTO((RL_VOID *params));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_FilteredStat (filtered.c)
*$ Abstract:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*$ Keywords:
*	PROFILE, SERIES
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	PRO_OBJECT	*Pro_FilteredStat(stat, filter, df1, df2)
*	PRO_OBJECT	*stat;
*	RL_FLT8		*filter;
*	RL_INT4		df1, df2;
*$ Inputs:
*	stat		pointer to the statistical series to resample.
*	filter[] 	array of filter coefficients.  The array should be of
*			size |df2-df1|+1.
*	df1, df2	index offset of each array endpoint from filter center.
*			Typically, df1<0 and df2>0.
*$ Outputs:
*	none
*$ Returns:
*	pointer to a new statistical series, or NULL on non-fatal error.
*$ Detailed_description:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*
*	The new series is, in effect, a convolution of the old series with a
*	fixed filter.  The filter coefficients are always adjusted so that the
*	sum of the weights on each new sample is unity.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.  The new object links to the original statistical
*	series.
*$ Examples:
*	Suppose stat is a statistical series.  Then this code creates a new
*	statistical series in which each sample is an average of the three
*	nearest samples in the original series.
*
*	PRO_OBJECT	*stat3;
*	RL_FLT8		filter[3] = {1., 1., 1.};
*
*	stat3 = Pro_FilteredStat(stat, filter, -1, 1);
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if stat is NULL or is not a stat series.
*	RL_MEMORY_ERROR		on memory allocation failure.
*$ 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_FilteredStat(stat, filter, df1, df2)
PRO_OBJECT	*stat;
RL_FLT8		*filter;
RL_INT4		df1, df2;
{
ZPRO_FILTERED	*params;
PRO_OBJECT	*new;
RL_INT4		i, k1, k2;
RL_FLT8		x1, dx;

    /* Verify object type */
    if (XPro_StatPtr(stat) == NULL) return NULL;

    /* Create new filtered statistical series structure */
    params = (ZPRO_FILTERED *) XRL_Malloc(sizeof(ZPRO_FILTERED));
    if (params == NULL) return NULL;

    params->fsize = df2 - df1 + 1;
    if (df2 < df1) params->fsize = df1 - df2 + 1;

    params->filter = (RL_FLT8 *) XRL_Malloc(params->fsize * sizeof(RL_FLT8));
    if (params->filter == NULL) {
	XRL_Free((RL_VOID *) params);
	return NULL;
    }

    /* Fill structure */
    params->stat = stat;

    params->df1 = df1;
    params->df2 = df2;
    for (i=0; i < params->fsize; i++)
	params->filter[i] = filter[i];

    /* Reverse coefficients if necessary */
    if (df2 < df1) {
	params->df1 = df2;
	params->df2 = df1;
	for (i=0; i < params->fsize; i++)
	    params->filter[i] = filter[params->fsize - i];
    }

    /* Create new weighted statistical series object */
    Pro_SeriesSampling(stat, &x1, NULL, &dx);
    Pro_SeriesIndices(stat, &k1, &k2);

    new = Pro_WeightedStat(stat, k1, k2, x1, dx,
                           ZPro_FilteredKRange, ZPro_FilteredWeight,
                           params, sizeof(ZPRO_FILTERED), TRUE,
                           ZPro_PrintFiltered);
    XRL_Free((RL_VOID *) params);

/* NOTE: It is not necessary to enslave old object because Pro_WeightedStat()
   does that for us!
    if (new != NULL && XPro_EnslaveObject(new,stat)) {
	Pro_FreeObject(new);
	return NULL;
    }
*/

    /* Transfer coordinate names to new object */
    Pro_RenameObject(new, 1, Pro_ObjectName(stat,1));
    Pro_RenameObject(new, 2, Pro_ObjectName(stat,2));

    return new;
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_FilteredKRange(knew, k1, k2, params)
*
* This internal routine returns the range of nonzero weights in the old series
* corresponding to a given index in the new series.
*
* Inputs:
*	knew		new index.
*	params		pointer to the ZPRO_FILTERED structure.
*
* Outputs:
*	k1, k2		range of indices in old series for which weight is
*			nonzero.
*******************************************************************************/

static void ZPro_FilteredKRange(knew, k1, k2, params)
RL_INT4		knew, *k1, *k2;
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;

    filtered = (ZPRO_FILTERED *) params;

    *k1 = knew + filtered->df1;
    *k2 = knew + filtered->df2;
}

/*
********************************************************************************
* ZPro_FilteredWeight(knew, k, params)
*
* This internal routine returns the (un-normalized) weight of an old sample in
* the new series.
*
* Inputs:
*	knew		index in new series.
*	k		index in old series.
*	params		pointer to the ZPRO_FILTERED structure.
*
* Return:		weight value.
*******************************************************************************/

static RL_FLT8	ZPro_FilteredWeight(knew, k, params)
RL_INT4		knew, k;
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;
RL_INT4		index;

    filtered = (ZPRO_FILTERED *) params;

    index = k - knew - filtered->df1;
    if (index < 0 || index >= filtered->fsize) return 0.;

    return filtered->filter[index];
}

/*
********************************************************************************
* ZPro_PrintFiltered(params)
*
* This internal routine prints information about a filtered stat series.
*
* Inputs:
*	params		pointer to the ZPRO_FILTERED data structure.
*******************************************************************************/

static void	ZPro_PrintFiltered(params)
RL_VOID		*params;
{
ZPRO_FILTERED	*filtered;
RL_INT4		i;

    filtered = (ZPRO_FILTERED *) params;

    /* Print filtered statistical series info... */
    printf("\nFiltered statistical series parameters...\n");

    for (i = filtered->df1; i <= filtered->df2; i++) {
	printf("filter[%d] = %#g\n", i, filtered->filter[i - filtered->df1]);
    }
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
*	FPro_FilteredStat (filtered.c)
*$ Abstract:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*$ Keywords:
*	PROFILE, SERIES
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	integer*4 function FPro_FilteredStat(stat, filter, df1, df2)
*	integer*4	stat, df1, df2
*	real*8		filter(*)
*$ Inputs:
*	stat		FORTRAN pointer to the statistical series to resample.
*	filter		array of filter coefficients.  The array should be of
*			size |df2-df1|+1.
*	df1, df2	index offset of each array endpoint from filter center.
*			Typically, df1<0 and df2>0.
*$ Outputs:
*	none
*$ Returns:
*	FORTRAN pointer to a new statistical series, or 0 on non-fatal error.
*$ Detailed_description:
*	This routine creates and returns a filtered statistical series, where
*	each new sample is a weighted sum of samples in another series.
*
*	The new series is, in effect, a convolution of the old series with a
*	fixed filter.  The filter coefficients are always adjusted so that the
*	sum of the weights on each new sample is unity.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.  The new object links to the original statistical
*	series.
*$ Examples:
*	Suppose stat is a statistical series.  Then this code creates a new
*	statistical series in which each sample is an average of the three
*	nearest samples in the original series.
*
*	integer*4	stat3
*	real*8		filter(3)/ 1.d0, 1.d0, 1.d0 /
*
*	stat3 = FPro_FilteredStat(stat, filter, -1, 1)
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if stat is NULL or is not a stat series.
*	RL_MEMORY_ERROR		on memory allocation failure.
*	FORTRAN_POINTER_ERROR	if stat 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_filteredstat) (stat, filter, df1, df2)
RL_INT4	*stat, *df1, *df2;
RL_FLT8	*filter;
{
RL_VOID *ptr;
RL_INT4 index;

    /* Look up statistical series pointer */
    ptr = FORT_GetPointer(*stat);
    if (ptr == NULL) return 0;

    /* Create filtered statistical series */
    ptr = (RL_VOID *) Pro_FilteredStat((PRO_OBJECT *) ptr, filter, *df1, *df2);
    if (ptr == NULL) return 0;

    /* Return new pointer */
    index = FORT_AddPointer(ptr);
    if (index == 0) Pro_FreeObject((PRO_OBJECT *) ptr);

    return index;
}

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