Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions lib/aggregate.dart
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,12 @@ export 'src/aggregate/checkpoint.dart';
export 'src/aggregate/laws.dart';
export 'src/aggregate/primitives.dart';
export 'src/onehz/clinical/hrv_time_aggregator.dart' hide hrvTimeOneSlice;


export 'src/onehz/clinical/load_trimp_aggregator.dart';
export 'src/onehz/clinical/cosinor_aggregator.dart';
export 'src/onehz/sleep/hypnogram_aggregator.dart';

export 'src/onehz/sleep/circadian_np_aggregator.dart';
export 'src/onehz/clinical/prsa_aggregator.dart';
export 'src/onehz/foundations/rr_sandwich.dart';



export 'src/onehz/clinical/readiness_lnrmssd_aggregator.dart';
export 'src/onehz/sleep/sri_aggregator.dart';
export 'src/onehz/workout/time_in_zone_aggregator.dart';
3 changes: 2 additions & 1 deletion lib/onehz.dart
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@ export 'src/onehz/clinical/hrv_freq.dart';
export 'src/onehz/clinical/prsa.dart';
export 'src/onehz/clinical/nocturnal.dart';
export 'src/onehz/clinical/illness_cusum.dart';
export 'src/onehz/clinical/readiness_lnrmssd.dart';
export 'src/onehz/clinical/readiness_lnrmssd.dart'
hide readinessLnRmssdFromSuffix;
export 'src/onehz/clinical/cosinor.dart';
export 'src/onehz/clinical/load_trimp.dart';
export 'src/onehz/clinical/stress_si.dart';
Expand Down
156 changes: 10 additions & 146 deletions lib/src/onehz/clinical/cosinor.dart
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,11 @@
// the daytime peak — callers choose the period and supply trough-anchored
// times; this engine just fits. R² is reported so a poor fit is visible.

import 'dart:math' as math;
import 'dart:typed_data';
import '../../aggregate/aggregator.dart';
import '../types.dart';
import '../util.dart';
import 'cosinor_aggregator.dart';

class CosinorFit {
final double mesor; // rhythm-adjusted mean
Expand Down Expand Up @@ -60,7 +62,6 @@ class CosinorFit {
/// reported.
const double cosinorAcrophaseAlpha = 0.05;


/// Minimum samples for a cosinor fit. The model has 3 free parameters
/// (M, β, γ); with n = 4 there is a single residual degree of freedom, so R²
/// is essentially an interpolation score (4 RANDOM points routinely fit at
Expand All @@ -86,154 +87,17 @@ Metric<CosinorFit> cosinor(
note: 'invalid cosinor period',
);
}
// Runs [cosinorAggregator] over the whole series in one slice. Mismatched
// lengths cannot be a [TimedValues] batch; they are the same "too few
// points" refusal as before.
if (tHours.length != y.length) {
return const Metric<CosinorFit>.absent(tier: Tier.high,
inputs_used: ['signal_timeseries'], note: 'too few points for cosinor');
}
final state = _CosinorState();
final w = 2 * math.pi / periodHours;
for (var i = 0; i < y.length; i++) {
state.addCentered(math.cos(w * tHours[i]), math.sin(w * tHours[i]), y[i]);
}
return _finishCosinor(state, periodHours: periodHours, minPoints: minPoints);
}

// Nelson et al. (1979), Methods for cosinor-rhythmometry.
// Centered least squares changes rounding. Relative 1e-12 gate is a
// numerical identifiability policy, not a physiological cutoff.
class _CosinorState {
var n = 0;
var meanC = 0.0;
var meanS = 0.0;
var meanY = 0.0;
var mCC = 0.0;
var mSS = 0.0;
var mYY = 0.0;
var mCS = 0.0;
var mCY = 0.0;
var mSY = 0.0;

void addCentered(double c, double s, double y) {
final nextN = n + 1;
final dc = c - meanC;
final ds = s - meanS;
final dy = y - meanY;
final nextMeanC = meanC + dc / nextN;
final nextMeanS = meanS + ds / nextN;
final nextMeanY = meanY + dy / nextN;
mCC += dc * (c - nextMeanC);
mSS += ds * (s - nextMeanS);
mYY += dy * (y - nextMeanY);
mCS += dc * (s - nextMeanS);
mCY += dc * (y - nextMeanY);
mSY += ds * (y - nextMeanY);
meanC = nextMeanC;
meanS = nextMeanS;
meanY = nextMeanY;
n = nextN;
}

}

Metric<CosinorFit> _finishCosinor(
_CosinorState state, {
required double periodHours,
required int minPoints,
}) {
const inputs = ['signal_timeseries'];
final n = state.n;
if (n < math.max(4, minPoints) || periodHours <= 0) {
return const Metric<CosinorFit>.absent(
tier: Tier.high,
inputs_used: inputs,
inputs_used: ['signal_timeseries'],
note: 'too few points for cosinor',
);
}
// A finite period or response can still overflow the design arithmetic.
// Reject invalid moments before NaN can bypass the identifiability gate.
if (![state.meanC, state.meanS, state.meanY, state.mCC, state.mSS,
state.mYY, state.mCS, state.mCY, state.mSY].every((v) => v.isFinite)) {
return const Metric<CosinorFit>.absent(
tier: Tier.high,
inputs_used: inputs,
note: 'non-finite cosinor design',
);
}
final determinant = state.mCC * state.mSS - state.mCS * state.mCS;
// The determinant is relative to the diagonal co-moments. A design needs
// at least 1e-12 of that scale after removing its means to identify both
// cosine coefficients; below that, numerical noise decides the fit.
if (determinant <= 1e-12 * state.mCC * state.mSS) {
return const Metric<CosinorFit>.absent(
tier: Tier.high,
inputs_used: inputs,
note: 'singular cosinor design (insufficient phase coverage)',
);
}
final beta =
(state.mCY * state.mSS - state.mSY * state.mCS) / determinant;
final gamma =
(state.mSY * state.mCC - state.mCY * state.mCS) / determinant;
final mesor = state.meanY - beta * state.meanC - gamma * state.meanS;
final amplitude = math.sqrt(beta * beta + gamma * gamma);
final acro = math.atan2(-gamma, beta);
final w = 2 * math.pi / periodHours;
var phaseHours = -acro / w;
phaseHours = phaseHours % periodHours;
if (phaseHours < 0) phaseHours += periodHours;

final ssTot = state.mYY;
final ssRes = state.mYY -
2 * beta * state.mCY -
2 * gamma * state.mSY +
beta * beta * state.mCC +
gamma * gamma * state.mSS +
2 * beta * gamma * state.mCS;
final double r2 = ssTot == 0 ? 0.0 : (1 - ssRes / ssTot).clamp(0, 1);
final double r2Adj = (1 - (1 - r2) * (n - 1) / (n - 3)).clamp(0, 1);
// Zero-amplitude F-test (Nelson 1979): the acrophase is only defined once
// the test rejects. A flat or rhythm-free series still yields some
// atan2(−γ, β) from residual noise or float rounding.
final dof = n - 3;
final double ampP;
if (ssTot <= 0) {
ampP = 1.0;
} else if (ssRes <= 0) {
ampP = 0.0;
} else {
final f = ((ssTot - ssRes) / 2) / (ssRes / dof);
ampP = _f2Survival(f, dof);
}
final phaseDefined = ampP < cosinorAcrophaseAlpha;

return Metric<CosinorFit>(
value: CosinorFit(
mesor: mesor,
amplitude: amplitude,
acrophaseRad: phaseDefined ? acro : null,
acrophaseHours: phaseDefined ? phaseHours : null,
periodHours: periodHours,
r2: r2,
r2Adj: r2Adj,
amplitudeP: ampP,
),
confidence: r2Adj.clamp(0.1, 0.95),
tier: Tier.high,
inputs_used: inputs,
note: phaseDefined
? 'single-component cosinor; anchor circadian phase on the nocturnal '
'trough, not the daytime peak'
: 'single-component cosinor; no resolvable rhythm '
'(zero-amplitude test not rejected), acrophase undefined',
);
}


/// Upper tail of F(2, d): P(F > f) = (1 + 2f/d)^(−d/2), exact for 2 numerator
/// degrees of freedom.
double _f2Survival(double f, int d) {
if (f.isNaN) return 1.0;
if (f.isInfinite) return 0.0;
if (f <= 0) return 1.0;
return math.pow(1 + 2 * f / d, -d / 2).toDouble();
return cosinorAggregator(periodHours: periodHours, minPoints: minPoints)
.run(TimedValues(
Float64List.fromList(tHours), Float64List.fromList(y)));
}
Loading
Loading