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 described in 4.1 Per-Point State and Update.

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 26 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 41 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: