/* weighted.c
********************************************************************************
* weighted.c -- Routines for re-weighted statistical series objects
*
* User routines
*	Pro_WeightedStat()	creates a re-weighted statistical series object.
*
* In a re-weighted statistical series, each sample is derived by taking a
* weighted sum of the samples in another statistical series.  The user must
* provide a function that calculates the (possibly un-normalized) weights and
* another function that returns the range of nonzero-weighted samples in the old
* series that comprise each new sample.
*
* Version 1.0: Original version.
*              Mark Showalter, PDS Rings Node, March 1998.
* Version 1.1: Support for Pro_StatRange() function added.
*              Mark Showalter, January 2000.
*******************************************************************************/
#include 
#include 
#include "profile.h"

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

typedef struct ZPRO_WEIGHTED_STRUCT {
    XPRO_CLASS	class;
    PRO_OBJECT	*oldstat;
    void	(*krange)  RL_PROTO((RL_INT4 knew, RL_INT4 *k1, RL_INT4 *k2,
		                     RL_VOID *params));
    RL_FLT8	(*weights) RL_PROTO((RL_INT4 knew, RL_INT4 k,
		                     RL_VOID *params));
    RL_VOID	*params;
    RL_SIZE	paramsize;
    RL_BOOL	isduped;
    RL_INT4	kmin, kmax;
    void	(*printmore) RL_PROTO((RL_VOID *params));
} ZPRO_WEIGHTED;

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

static XPRO_CLASS weighted_class = {XPRO_WEIGHTED_CLASS, "weighted", NULL};

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

static RL_FLT8	ZPro_WeightedValue RL_PROTO((RL_VOID *pointer, RL_INT4 k,
		                             RL_INT4 *flag));
static RL_FLT8	ZPro_WeightedCovar RL_PROTO((RL_VOID *pointer, RL_INT4 k1,
		                             RL_INT4 k2));
static void	ZPro_WeightedRange RL_PROTO((RL_VOID *pointer, RL_INT4 k,
		                             RL_INT4 *k1, RL_INT4 *k2));
static void 	ZPro_FreeWeighted  RL_PROTO((RL_VOID *pointer));
static void	ZPro_PrintWeighted RL_PROTO((RL_VOID *pointer));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_WeightedStat (weighted.c)
*$ Abstract:
*	This routine creates a new statistical series in which each sample is a
*	weighted sum of samples from another statistical series.
*$ Keywords:
*	PROFILE
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	PRO_OBJECT	*Pro_WeightedStat(stat, k1, k2, x1, dx, krange, weights,
*		                 params, paramsize, dupe, printmore)
*	PRO_OBJECT	*stat;
*	RL_INT4		k1, k2;
*	RL_FLT8		x1, dx;
*	void		(*krange) (RL_INT4 knew, RL_INT4 *k1, RL_INT4 *k2,
*			           RL_VOID *params);
*	RL_FLT8		(*weights) (RL_INT4 knew, RL_INT4 k, RL_VOID *params);
*	RL_VOID		*params;
*	RL_SIZE		paramsize;
*	RL_BOOL		dupe;
*	void		(*printmore) (RL_VOID *params);
*$ Inputs:
*	stat		pointer to the statistical series to re-weight.
*	k1, k2		new range of indices.
*	x1, dx		new x coordinate starting value and interval.
*	krange		function that returns the range of indices on the old
*			series corresponding to nonzero weights for a given new
*			index.
*	weights		function that returns the weights on each old sample
*			for each new sample.
*	params		parameter set to pass to krange and weight functions.
*	paramsize	number of bytes in parameter set.
*	dupe		TRUE to duplicate parameters; FALSE just to save
*			pointer.
*	printmore	optional subroutine to print more info about weighted
*			stat series; NULL to ignore.
*$ Outputs:
*	none
*$ Returns:
*	pointer to a new statistical series object; NULL on non-fatal error.
*$ Detailed_description:
*	This routine creates a new statistical series in which each sample is a
*	weighted sum of samples from another statistical series.  The user must
*	provide subroutines that return the set of samples and (possibly
*	unnormalized) weights to apply to each old sample to generate any given
*	new sample.
*
*	This is the "workhorse" routine used by Pro_Filtered() and
*	Pro_Resampled().  Note that there is no FORTRAN-callable equivalent to
*	this routine.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.  The new object links to the original statistical
*	series.
*$ Examples:
*	Suppose we have a statistical series stat, and we wish to create a new
*	series in which each sample is the average of the nearest three samples
*	in the original series.  (This is actually much easier to do using
*	Pro_FilteredStat).
*
*	We need these two functions:
*
*	void krange(RL_INT4 knew, RL_INT4 *k1, RL_INT4 *k2, RL_VOID *params)
*	{
*		k1 = knew - 1;
*		k2 = knew + 1;
*	}
*
*	RL_FLT8	weights(RL_INT4 knew, RL_INT4 k, RL_VOID *params)
*	{
*		return 1.;
*	}
*
*	Then this snippet of code creates and returns stat3, the new series:
*
*	PRO_OBJECT	*stat, *stat3;
*	RL_INT4		k1, k2;
*	RL_FLT8		x1, dx;

*	Pro_SeriesIndices(stat, &k1, &k2);
*	Pro_SeriesSampling(stat, &x1, NULL, &dx);
*	stat3 = Pro_WeightedStat(stat, k1, k2, x1, dx, krange, weights,
*			NULL, 0, FALSE, NULL);
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	RL_MEMORY_ERROR		on memory allocation failure.
*	PRO_CLASS_ERROR 	if the object is NULL or is not a stat series.
*$ 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_WeightedStat(stat, k1, k2, x1, dx, krange, weights,
		                 params, paramsize, dupe, printmore)
PRO_OBJECT	*stat;
RL_INT4		k1, k2;
RL_FLT8		x1, dx;
void		(*krange)  RL_PROTO((RL_INT4 knew, RL_INT4 *k1, RL_INT4 *k2,
		                     RL_VOID *params));
RL_FLT8		(*weights) RL_PROTO((RL_INT4 knew, RL_INT4 k,
		                     RL_VOID *params));
RL_VOID		*params;
RL_SIZE		paramsize;
RL_BOOL		dupe;
void		(*printmore) RL_PROTO((RL_VOID *params));
{
ZPRO_WEIGHTED	*weighted;
PRO_OBJECT	*new;

    /* Make sure it's a statistical series */
    if (XPro_StatPtr(stat) == NULL) return NULL;

    /* Allocate new weighted statistical series structure */
    weighted = (ZPRO_WEIGHTED *) XRL_Malloc(sizeof(ZPRO_WEIGHTED));
    if (weighted == NULL) return NULL;

    /* Initialize structure */
    weighted->class     = weighted_class;
    weighted->oldstat   = stat;
    weighted->krange    = krange;
    weighted->weights   = weights;
    weighted->printmore = printmore;
    weighted->kmin      = k1;
    weighted->kmax      = k2;

    /* Copy parameters if necessary */
    if (dupe && paramsize > 0 && params != NULL) {
	weighted->params = XRL_Malloc(paramsize);
	if (weighted->params == NULL) {
	    XRL_Free((RL_VOID *) weighted);
	    return NULL;
	}
	memcpy(weighted->params, params, paramsize);
	weighted->paramsize = paramsize;
	weighted->isduped   = TRUE;
    }
    else {
	weighted->params    = params;
	weighted->paramsize = paramsize;
	weighted->isduped   = FALSE;
    }

    /* Create new statistical series object */
    new = XPro_MakeStat(k1, k2, x1, dx,
                        ZPro_WeightedValue,
			ZPro_WeightedCovar, ZPro_WeightedRange,
                        ZPro_FreeWeighted, ZPro_PrintWeighted,
                        (RL_VOID *) weighted);

    /* Enslave old statistical series and check for errors */
    if (new != NULL && XPro_EnslaveObject(new, stat)) {
	Pro_FreeObject(new);
	new = NULL;
    }

    /* Transfer Y-coordinate name to new object */
    Pro_RenameObject(new, 2, Pro_ObjectName(stat,2));

    return new;
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_WeightedValue(pointer, k, flag)
*
* This internal routine returns the value of a weighted statistical series.
*
* Inputs:
*	pointer		pointer to the ZPRO_WEIGHTED structure.
*	k		index at which to evaluate series.
*
* Outputs:
*	*flag		0 if value is valid.
*******************************************************************************/

static RL_FLT8	ZPro_WeightedValue(pointer, k, flag)
RL_VOID		*pointer;
RL_INT4		k, *flag;
{
ZPRO_WEIGHTED	*stat;
RL_INT4		i, i1, i2, oldflag, newflag;
RL_FLT8		sum0, sum1, value, weight;

    stat = (ZPRO_WEIGHTED *) pointer;
    (stat->krange) (k, &i1, &i2, stat->params);

    sum0 = 0.;
    sum1 = 0.;
    newflag = PRO_MISSING_FLAG;	/* use only if every value is missing */
    for (i=i1; i<=i2; i++) {
	value = Pro_SeriesValue(stat->oldstat, i, &oldflag);

 	if (oldflag == 0) {
	    weight = (stat->weights)(k, i, stat->params);
	    sum0 = sum0 + weight;
	    sum1 = sum1 + weight * value;
	}
	else if (oldflag == PRO_INVALID_FLAG) {
	    newflag = PRO_INVALID_FLAG; /* change from missing to invalid */
	}
    }

    if (sum0 == 0.) {
	*flag = newflag;
	return 0.;
    }

    *flag = 0;
    return sum1/sum0;
}

/*
********************************************************************************
* ZPro_WeightedCovar(pointer, k1, k2)
*
* This internal routine evaluates the covariance of a weighted statistical
* series.
*
* Inputs:
*	pointer		pointer to the ZPRO_WEIGHTED structure.
*	k1, k2		pair of new indices at which to evaluate covariance.
*
* Returns:		value of covariance between k1 and k2;
*******************************************************************************/

static RL_FLT8	ZPro_WeightedCovar(pointer, k1, k2)
RL_VOID		*pointer;
RL_INT4		k1, k2;
{
ZPRO_WEIGHTED	*stat;
RL_INT4		i, i1, i2, j, j1, j2, jmin, jmax, flag;
RL_FLT8		sum0, sum1, wsumi, wsumj, weight;

    stat = (ZPRO_WEIGHTED *) pointer;

    /* Find ranges of associated samples in old series */
    (stat->krange) (k1, &i1, &i2, stat->params);
    (stat->krange) (k2, &j1, &j2, stat->params);

/*******************************************************************************
** The old code...
**
**  sum0 = 0.;
**  sum1 = 0.;
**  for (i=i1; i<=i2; i++) {
**	(void) Pro_SeriesValue(stat->oldstat, i, &flag1);
**
**  for (j=j1; j<=j2; j++) {
**	(void) Pro_SeriesValue(stat->oldstat, j, &flag2);
**
**	if (flag1 == 0 && flag2 == 0) {
**	    weight = (stat->weights)(k1, i, stat->params) *
**	             (stat->weights)(k2, j, stat->params);
**	    sum0 = sum0 + weight;
**	    sum1 = sum1 + weight * Pro_StatCovar(stat->oldstat, i, j);
**	}
**  }}
**
** The trouble with the code above is that, as the filter width increases, the
** computational effort goes as the width squared.  It fails to take advantage
** of the fact that the covariance of samples i and j in the original series is
** likely to be zero except for i close to j.  The revised code shown below
** takes advantage of this fact and so runs much faster.
*******************************************************************************/

/* Sum the products of weights for each pair of old samples */
    wsumi = 0.;
    for (i=i1; i<=i2; i++) {
        (void) Pro_SeriesValue(stat->oldstat, i, &flag);
	if (flag == 0) wsumi += (stat->weights)(k1, i, stat->params);
    }

    if (k1 == k2) {
	wsumj = wsumi;
    } else {
	wsumj = 0.;
	for (j=j1; j<=j2; j++) {
            (void) Pro_SeriesValue(stat->oldstat, j, &flag);
	    if (flag == 0) wsumj += (stat->weights)(k2, j, stat->params);
	}
    }

    sum0 = wsumi * wsumj;

/* Sum the products of weights times the covariance of each pair.
 * For this sum, we can skip pairs whose covariance is zero.
 */
    sum1 = 0.;
    for (i=i1; i<=i2; i++) {
	(void) Pro_SeriesValue(stat->oldstat, i, &flag);
	if (flag != 0) continue;

	Pro_StatRange(stat->oldstat, i, &jmin, &jmax);
	if (jmax < j1) continue;
	if (jmin > j2) continue;

	if (jmin < j1) jmin = j1;
	if (jmax > j2) jmax = j2;

	weight = (stat->weights)(k1, i, stat->params);
	for (j=jmin; j<=jmax; j++) {
	    (void) Pro_SeriesValue(stat->oldstat, j, &flag);
	    if (flag != 0) continue;

	    weight *= (stat->weights)(k2, j, stat->params);
	    sum1 += weight * Pro_StatCovar(stat->oldstat, i, j);
	}
    }

    if (sum0 == 0.) return 0.;

    return sum1/sum0;
}

/*
********************************************************************************
* ZPro_WeightedRange(pointer, k, k1, k2)
*
* This internal routine evaluates the range of samples that might have non-zero
* covariances with a given sample in a weighted statistical series.
*
* Inputs:
*	pointer		pointer to the ZPRO_WEIGHTED structure.
*	k		index at which to evaluate correlation range.
*
* Outputs:
*	k1, k2		range of indices within which samples are correlated.
*
* Returns:		none.
*******************************************************************************/

static void	ZPro_WeightedRange(pointer, k, k1, k2)
RL_VOID		*pointer;
RL_INT4		k, *k1, *k2;
{
ZPRO_WEIGHTED	*stat;
RL_INT4		i1, i2, imin, imax, dummy, ktest;

    stat = (ZPRO_WEIGHTED *) pointer;

    /* Find range of old samples that are used in this sample */
    (stat->krange) (k, &i1, &i2, stat->params);

    /* Find range of old samples that might correlate with this range */
    Pro_StatRange(stat->oldstat, i1, &imin,  &dummy);
    Pro_StatRange(stat->oldstat, i2, &dummy, &imax);

    /* Search down for first uncorrelated new sample */
    for (ktest = k; ktest >= stat->kmin; ktest--) {
	(stat->krange) (ktest, &i1, &i2, stat->params);
	if (i2 < imin || i1 > imax) break;
    }
    *k1 = ktest + 1;

    /* Search up for first uncorrelated sample */
    for (ktest = k; ktest <= stat->kmax; ktest++) {
	(stat->krange) (ktest, &i1, &i2, stat->params);
	if (i2 < imin || i1 > imax) break;
    }
    *k2 = ktest - 1;
}

/*
********************************************************************************
* ZPro_FreeWeighted(pointer)
*
* This internal routine deallocates the memory used by a weighted statistical
* series.
*
* Input:
*	pointer		pointer to the ZPRO_WEIGHTED structure.
*******************************************************************************/

static void	ZPro_FreeWeighted(pointer)
RL_VOID		*pointer;
{
ZPRO_WEIGHTED	*stat;

    stat = (ZPRO_WEIGHTED *) pointer;

    if (stat->isduped) XRL_Free(stat->params);
    XRL_Free(pointer);
}

/*
********************************************************************************
* ZPro_PrintWeighted(pointer)
*
* This internal routine prints information about a weighted statistical series.
*
* Input:
*	pointer		pointer to the ZPRO_FIXEDSTAT structure.
*******************************************************************************/

static void	ZPro_PrintWeighted(pointer)
RL_VOID		*pointer;
{
ZPRO_WEIGHTED	*weighted;

    weighted = (ZPRO_WEIGHTED *) pointer;

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

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

    /* Print object info... */
    printf("\nWeighted statistical series parameters...\n");
    printf("stat series = "); XPro_PrintInfo(weighted->oldstat);
    printf("  paramsize = %d\n", weighted->paramsize);
    printf(" duplicated = %s\n", (weighted->isduped ? "true" : "false"));

    /* Print optional additional info... */
    if (weighted->printmore != NULL)
	(*weighted->printmore) (weighted->params);

    return;
}

/*******************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
* Not provided, because generating FORTRAN interfaces to the "krange" and
* "weights" functions is just too complicated.
********************************************************************************
*/