/* 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 Ring-Moon Systems Node, March 1998.
* Version 1.1: Support for Pro_StatRange() function added.
* Mark Showalter, January 2000.
*******************************************************************************/
#include <stdio.h>
#include <string.h>
#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.
********************************************************************************
*/