statistics_moments.h File Reference

Weighted centered-moment kernels for the field-statistics pipeline. More...

#include <petscsys.h>

Include dependency graph for statistics_moments.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Data Structures
struct	PicurvMomentState
	Weighted centered state for one scalar quantity at one point. More...
struct	PicurvCoMomentState
	Weighted centered co-moment state for one ordered pair of quantities. More...

Functions
void	PicurvMomentStateReset (PicurvMomentState *state)
	Resets a scalar moment accumulator to the empty state.
void	PicurvCoMomentStateReset (PicurvCoMomentState *state)
	Resets a co-moment accumulator to the empty state.
PetscErrorCode	PicurvMomentStateUpdate (PicurvMomentState *state, PetscReal value, PetscReal weight)
	Applies one weighted sample to a scalar moment accumulator.
PetscErrorCode	PicurvCoMomentStateUpdate (PicurvCoMomentState *state, PetscReal value_x, PetscReal value_y, PetscReal weight)
	Applies one weighted paired sample to a co-moment accumulator.
PetscErrorCode	PicurvMomentStateMerge (PicurvMomentState *result, const PicurvMomentState *a, const PicurvMomentState *b)
	Merges two independently accumulated scalar states.
PetscErrorCode	PicurvCoMomentStateMerge (PicurvCoMomentState *result, const PicurvCoMomentState *a, const PicurvCoMomentState *b)
	Merges two independently accumulated co-moment states.
PetscReal	PicurvMomentStateVariance (const PicurvMomentState *state)
	Returns the weighted variance M2/W, or zero when no weight accumulated.
PetscReal	PicurvCoMomentStateCovariance (const PicurvCoMomentState *state)
	Returns the weighted covariance C/W, or zero when no weight accumulated.
PetscReal	PicurvMomentStateEffectiveCount (const PicurvMomentState *state)
	Returns Kish effective sample size W^2/W2, or zero when no weight accumulated.

Detailed Description

Weighted centered-moment kernels for the field-statistics pipeline.

These are pure numerical kernels: no configuration, no PETSc vectors, and no knowledge of windows, fields, or scheduling. They implement the accumulator contract fixed in Field Statistics Phase 2 Implementation Specification section 5, which in turn follows Field Statistics Pipeline Specification section 7.

Centered state is used rather than raw sums because it retains everything the first and second moments need while staying numerically stable for high-mean/low-fluctuation signals. Raw sums are recoverable as S_x = W*mu and Q_xx = M2 + W*mu^2.

Definition in file statistics_moments.h.

---

Data Structure Documentation

PicurvMomentState

struct PicurvMomentState

Weighted centered state for one scalar quantity at one point.

weight_sq is retained because physical-time weighting produces unequal weights, so an effective sample size cannot be inferred from count alone.

Definition at line 28 of file statistics_moments.h.

Data Fields
PetscReal	count	Number of accepted samples.
PetscReal	weight	Total weight W.
PetscReal	weight_sq	Sum of squared weights W2.
PetscReal	mean	Weighted mean mu.
PetscReal	m2	Centered second-moment sum M2.

PicurvCoMomentState

struct PicurvCoMomentState

Weighted centered co-moment state for one ordered pair of quantities.

The two means are tracked inside the pair rather than referenced from two PicurvMomentState values, so a co-moment update is self-contained and cannot be corrupted by update ordering between the members.

Definition at line 43 of file statistics_moments.h.

Data Fields
PetscReal	count	Number of accepted samples.
PetscReal	weight	Total weight W.
PetscReal	weight_sq	Sum of squared weights W2.
PetscReal	mean_x	Weighted mean of the first member.
PetscReal	mean_y	Weighted mean of the second member.
PetscReal	cm	Centered co-moment sum C.

Function Documentation

PicurvMomentStateReset()

void PicurvMomentStateReset

(

PicurvMomentState *

state

)

Resets a scalar moment accumulator to the empty state.

Parameters

[out]

state

Accumulator to clear; ignored when NULL.

Resets a scalar moment accumulator to the empty state.

See also

PicurvMomentStateReset()

Definition at line 15 of file statistics_moments.c.

{
    if (state == NULL) return;
    state->count = 0.0;
    state->weight = 0.0;
    state->weight_sq = 0.0;
    state->mean = 0.0;
    state->m2 = 0.0;
}

Here is the caller graph for this function:

PicurvCoMomentStateReset()

void PicurvCoMomentStateReset

(

PicurvCoMomentState *

state

)

Resets a co-moment accumulator to the empty state.

Parameters

[out]

state

Accumulator to clear; ignored when NULL.

Resets a co-moment accumulator to the empty state.

See also

PicurvCoMomentStateReset()

Definition at line 29 of file statistics_moments.c.

{
    if (state == NULL) return;
    state->count = 0.0;
    state->weight = 0.0;
    state->weight_sq = 0.0;
    state->mean_x = 0.0;
    state->mean_y = 0.0;
    state->cm = 0.0;
}

Here is the caller graph for this function:

PicurvMomentStateUpdate()

PetscErrorCode PicurvMomentStateUpdate	(	PicurvMomentState *	state,
		PetscReal	value,
		PetscReal	weight
	)

Applies one weighted sample to a scalar moment accumulator.

Implements the stable weighted update ‘W’ = W+w,d = x-mu,mu' = mu + (w/W')d,M2' = M2 + w*d*(x-mu')`.

Parameters

[in,out]	state	Accumulator to update.
[in]	value	Sample value x.
[in]	weight	Sample weight w; must be strictly positive.

Returns

Zero on success, or PETSC_ERR_ARG_OUTOFRANGE for a non-positive weight.

Applies one weighted sample to a scalar moment accumulator.

See also

PicurvMomentStateUpdate()

Definition at line 44 of file statistics_moments.c.

{
    PetscReal new_weight = 0.0;
    PetscReal delta = 0.0;
 
    PetscFunctionBeginUser;
    PetscCheck(state != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Moment accumulator is required.");
    PetscCheck(weight > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Sample weight must be strictly positive, got %g.", (double)weight);
 
    new_weight = state->weight + weight;
    delta = value - state->mean;
    state->mean += (weight / new_weight) * delta;
    /* Uses the post-update mean deliberately: this is the pairing that keeps the
     * centered sum stable when the mean is large relative to the fluctuation. */
    state->m2 += weight * delta * (value - state->mean);
    state->weight = new_weight;
    state->weight_sq += weight * weight;
    state->count += 1.0;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvCoMomentStateUpdate()

PetscErrorCode PicurvCoMomentStateUpdate	(	PicurvCoMomentState *	state,
		PetscReal	value_x,
		PetscReal	value_y,
		PetscReal	weight
	)

Applies one weighted paired sample to a co-moment accumulator.

Uses the pre-update mean of x and the post-update mean of y, which is the form that reduces exactly to the M2 update when x and y are the same signal.

Parameters

[in,out]	state	Accumulator to update.
[in]	value_x	Sample value of the first member.
[in]	value_y	Sample value of the second member.
[in]	weight	Sample weight w; must be strictly positive.

Returns

Zero on success, or PETSC_ERR_ARG_OUTOFRANGE for a non-positive weight.

Applies one weighted paired sample to a co-moment accumulator.

See also

PicurvCoMomentStateUpdate()

Definition at line 71 of file statistics_moments.c.

{
    PetscReal new_weight = 0.0;
    PetscReal delta_x = 0.0;
 
    PetscFunctionBeginUser;
    PetscCheck(state != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Co-moment accumulator is required.");
    PetscCheck(weight > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
               "Sample weight must be strictly positive, got %g.", (double)weight);
 
    new_weight = state->weight + weight;
    delta_x = value_x - state->mean_x;
    state->mean_x += (weight / new_weight) * delta_x;
    state->mean_y += (weight / new_weight) * (value_y - state->mean_y);
    /* Pre-update mean of x against post-update mean of y. With x == y this
     * reduces exactly to the M2 update, which the unit tests assert. */
    state->cm += weight * delta_x * (value_y - state->mean_y);
    state->weight = new_weight;
    state->weight_sq += weight * weight;
    state->count += 1.0;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvMomentStateMerge()

PetscErrorCode PicurvMomentStateMerge	(	PicurvMomentState *	result,
		const PicurvMomentState *	a,
		const PicurvMomentState *	b
	)

Merges two independently accumulated scalar states.

Implements the stable weighted parallel combination d = mu_b - mu_a, W = W_a + W_b, mu = mu_a + (W_b/W)d, M2 = M2_a + M2_b + (W_a*W_b/W)d^2. Merging an empty state is a no-op, so partitions that received no samples are safe to combine.

Parameters

[out]	result	Merged state; may alias a or b.
[in]	a	First partition.
[in]	b	Second partition.

Returns

Zero on success, or PETSC_ERR_ARG_NULL when any argument is NULL.

Merges two independently accumulated scalar states.

See also

PicurvMomentStateMerge()

Definition at line 100 of file statistics_moments.c.

{
    PetscReal total_weight = 0.0;
    PetscReal delta = 0.0;
    PicurvMomentState merged;
 
    PetscFunctionBeginUser;
    PetscCheck(result != NULL && a != NULL && b != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Merge requires a destination and two source accumulators.");
 
    if (b->weight == 0.0) { *result = *a; PetscFunctionReturn(0); }
    if (a->weight == 0.0) { *result = *b; PetscFunctionReturn(0); }
 
    total_weight = a->weight + b->weight;
    delta = b->mean - a->mean;
    merged.count = a->count + b->count;
    merged.weight = total_weight;
    merged.weight_sq = a->weight_sq + b->weight_sq;
    merged.mean = a->mean + (b->weight / total_weight) * delta;
    merged.m2 = a->m2 + b->m2 + (a->weight * b->weight / total_weight) * delta * delta;
    *result = merged;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvCoMomentStateMerge()

PetscErrorCode PicurvCoMomentStateMerge	(	PicurvCoMomentState *	result,
		const PicurvCoMomentState *	a,
		const PicurvCoMomentState *	b
	)

Merges two independently accumulated co-moment states.

Parameters

[out]	result	Merged state; may alias a or b.
[in]	a	First partition.
[in]	b	Second partition.

Returns

Zero on success, or PETSC_ERR_ARG_NULL when any argument is NULL.

Merges two independently accumulated co-moment states.

See also

PicurvCoMomentStateMerge()

Definition at line 129 of file statistics_moments.c.

{
    PetscReal total_weight = 0.0;
    PetscReal delta_x = 0.0;
    PetscReal delta_y = 0.0;
    PicurvCoMomentState merged;
 
    PetscFunctionBeginUser;
    PetscCheck(result != NULL && a != NULL && b != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Merge requires a destination and two source accumulators.");
 
    if (b->weight == 0.0) { *result = *a; PetscFunctionReturn(0); }
    if (a->weight == 0.0) { *result = *b; PetscFunctionReturn(0); }
 
    total_weight = a->weight + b->weight;
    delta_x = b->mean_x - a->mean_x;
    delta_y = b->mean_y - a->mean_y;
    merged.count = a->count + b->count;
    merged.weight = total_weight;
    merged.weight_sq = a->weight_sq + b->weight_sq;
    merged.mean_x = a->mean_x + (b->weight / total_weight) * delta_x;
    merged.mean_y = a->mean_y + (b->weight / total_weight) * delta_y;
    merged.cm = a->cm + b->cm + (a->weight * b->weight / total_weight) * delta_x * delta_y;
    *result = merged;
    PetscFunctionReturn(0);
}

Here is the caller graph for this function:

PicurvMomentStateVariance()

PetscReal PicurvMomentStateVariance

(

const PicurvMomentState *

state

)

Returns the weighted variance M2/W, or zero when no weight accumulated.

This is the population (weight-normalized) variance the postprocessing contract uses for R_ii. It is never negative for a state built only through PicurvMomentStateUpdate, but callers taking a square root must still clamp, because cancellation can drive M2 slightly negative for degenerate inputs.

Parameters

[in]

state

Accumulator to read.

Returns

Weighted variance, or zero for a NULL or unsampled accumulator.

Returns the weighted variance M2/W, or zero when no weight accumulated.

See also

PicurvMomentStateVariance()

Definition at line 161 of file statistics_moments.c.

{
    if (state == NULL || state->weight == 0.0) return 0.0;
    return state->m2 / state->weight;
}

Here is the caller graph for this function:

PicurvCoMomentStateCovariance()

PetscReal PicurvCoMomentStateCovariance

(

const PicurvCoMomentState *

state

)

Returns the weighted covariance C/W, or zero when no weight accumulated.

Parameters

[in]

state

Accumulator to read.

Returns

Weighted covariance, or zero for a NULL or unsampled accumulator.

Returns the weighted covariance C/W, or zero when no weight accumulated.

See also

PicurvCoMomentStateCovariance()

Definition at line 171 of file statistics_moments.c.

{
    if (state == NULL || state->weight == 0.0) return 0.0;
    return state->cm / state->weight;
}

Here is the caller graph for this function:

PicurvMomentStateEffectiveCount()

PetscReal PicurvMomentStateEffectiveCount

(

const PicurvMomentState *

state

)

Returns Kish effective sample size W^2/W2, or zero when no weight accumulated.

Equals the sample count exactly under equal weights, and degrades toward one as the weight distribution becomes dominated by a single sample.

Parameters

[in]

state

Accumulator to read.

Returns

Effective sample size, or zero for a NULL or unsampled accumulator.

Returns Kish effective sample size W^2/W2, or zero when no weight accumulated.

See also

PicurvMomentStateEffectiveCount()

Definition at line 181 of file statistics_moments.c.

{
    if (state == NULL || state->weight_sq == 0.0) return 0.0;
    return (state->weight * state->weight) / state->weight_sq;
}

Here is the caller graph for this function: