diff --git a/src/lib/analytics/mood-prognosis/__tests__/features.test.ts b/src/lib/analytics/mood-prognosis/__tests__/features.test.ts new file mode 100644 index 000000000..19379217c --- /dev/null +++ b/src/lib/analytics/mood-prognosis/__tests__/features.test.ts @@ -0,0 +1,281 @@ +/** + * The matrix: what becomes a column, what never does, and what happens to a + * day that is missing something. + * + * The cases that matter are the absences. A day with no context row must not + * contribute a zero to a work-status column, because zero there means "worked + * a regular day and not overtime" and the truth is "nobody said". That + * distinction is invisible in the fitted numbers and would quietly teach the + * model that every unlogged day was a day off. + */ +import { describe, expect, it } from "vitest"; + +import { + MAX_MODEL_FEATURES, + MIN_FEATURE_COVERAGE, + buildFeatureMatrix, + extractDayFeatures, + screenFeatures, + standardiseDay, + type PrognosisContextInput, + type PrognosisDayInput, +} from "../features"; +import { MOOD_DIMENSIONS } from "@/lib/mood/dimensions"; + +const EMPTY_LINKED = { + sleepAsleep: null, + steps: null, + activeEnergy: null, + restingHeartRate: null, + heartRateVariability: null, +}; + +function context( + overrides: Partial = {}, +): PrognosisContextInput { + return { + workStatus: null, + workMinutes: null, + overtimeMinutes: null, + workLoad: null, + workSatisfaction: null, + contactCircles: [], + contactForm: null, + contactExtent: null, + contactQuality: null, + contactSupport: null, + leisureCategories: [], + leisureMinutes: null, + leisureJoy: null, + leisureRecovery: null, + eventType: null, + eventValence: null, + ...overrides, + }; +} + +function day( + index: number, + overrides: Partial = {}, +): PrognosisDayInput { + const date = new Date(Date.UTC(2026, 0, 1 + index)); + return { + day: date.toISOString().slice(0, 10), + a1: 5, + dimensions: {}, + context: null, + linked: { ...EMPTY_LINKED }, + ...overrides, + }; +} + +describe("extractDayFeatures", () => { + it("takes A2 through the dimension table's inverse flag, not inline", () => { + const a2 = MOOD_DIMENSIONS.find((d) => d.key === "a2"); + expect(a2?.inverse, "A2 is the inverse dimension this case is about").toBe( + true, + ); + const features = extractDayFeatures( + day(0, { dimensions: { a2: 8, a3: 8 } }), + ); + // Stress 8 is a heavy day, so on an up-is-better axis it reads as 2. + expect(features.get("dimension:a2")).toBe(2); + // A3 is not inverse and is carried literally. + expect(features.get("dimension:a3")).toBe(8); + }); + + it("never makes A1 an input", () => { + const features = extractDayFeatures( + day(0, { a1: 7, dimensions: { a1: 7, a2: 3 } }), + ); + expect(features.has("dimension:a1")).toBe(false); + }); + + it("a day with no context contributes no context column at all", () => { + const features = extractDayFeatures(day(0, { context: null })); + expect( + [...features.keys()].filter((k) => k.startsWith("context:")), + ).toEqual([]); + }); + + it("an unanswered field contributes no column; an answered one contributes zeros and a one", () => { + const answered = extractDayFeatures( + day(0, { context: context({ workStatus: "overtime" }) }), + ); + expect(answered.get("context:workStatus=overtime")).toBe(1); + expect(answered.get("context:workStatus=off")).toBe(0); + + const unanswered = extractDayFeatures( + day(0, { context: context({ workLoad: 4 }) }), + ); + expect(unanswered.has("context:workStatus=overtime")).toBe(false); + expect(unanswered.has("context:workStatus=off")).toBe(false); + expect(unanswered.get("context:workLoad")).toBe(4); + }); + + it("a recorded zero is a value, and absence is absence", () => { + const zero = extractDayFeatures( + day(0, { linked: { ...EMPTY_LINKED, steps: 0 } }), + ); + expect(zero.get("linked:steps")).toBe(0); + const absent = extractDayFeatures(day(0)); + expect(absent.has("linked:steps")).toBe(false); + }); + + it("an empty multi-select says nothing rather than saying no to everything", () => { + const features = extractDayFeatures( + day(0, { context: context({ contactCircles: [] }) }), + ); + expect(features.has("context:contactCircles=partner")).toBe(false); + }); +}); + +describe("screening", () => { + /** A2 drives A1 exactly; everything else is constant or noise. */ + function drivenDays(count: number): PrognosisDayInput[] { + return Array.from({ length: count }, (_, i) => { + const stress = i % 11; + return day(i, { + // Up-is-better A2 is `10 - stress`, and the target follows it. + a1: 10 - stress, + dimensions: { a2: stress }, + }); + }); + } + + it("keeps a feature that tracks the target", () => { + const screened = screenFeatures(drivenDays(40)); + const a2 = screened.find((s) => s.key === "dimension:a2"); + expect(a2?.selected, "a perfectly tracking feature was screened out").toBe( + true, + ); + expect(a2?.pairs).toBe(40); + }); + + it("refuses a feature recorded on too few of the days, before testing it", () => { + const days = drivenDays(40); + // Recorded on a quarter of the days — under the coverage floor whatever + // it correlates with. + const covered = Math.floor(40 * MIN_FEATURE_COVERAGE) - 5; + for (let i = 0; i < covered; i++) { + days[i] = { ...days[i], linked: { ...EMPTY_LINKED, steps: 1000 + i } }; + } + const screened = screenFeatures(days); + const steps = screened.find((s) => s.key === "linked:steps"); + expect(steps?.selected).toBe(false); + expect(steps?.reason).toBe("coverage"); + }); + + it("refuses a one-hot seen fewer times than the occurrence floor", () => { + const days = drivenDays(40).map((d, i) => ({ + ...d, + // Everybody worked; two days were overtime. The column is present on + // every day, so coverage passes and the occurrence floor is what + // refuses it. + context: context({ workStatus: i < 2 ? "overtime" : "regular" }), + })); + const screened = screenFeatures(days); + const overtime = screened.find( + (s) => s.key === "context:workStatus=overtime", + ); + expect(overtime?.selected).toBe(false); + expect(overtime?.reason).toBe("occurrences"); + }); + + it("refuses pure noise — the whole reason the FDR layer is here", () => { + // A deterministic pseudo-random series, so the case cannot flake. + let seed = 42; + const random = () => { + seed = (seed * 1103515245 + 12345) % 2147483648; + return seed / 2147483648; + }; + const days = Array.from({ length: 60 }, (_, i) => + day(i, { + a1: Math.round(random() * 10), + dimensions: { a2: Math.round(random() * 10) }, + linked: { ...EMPTY_LINKED, steps: Math.round(random() * 10000) }, + }), + ); + const selected = screenFeatures(days).filter((s) => s.selected); + expect(selected.map((s) => s.key)).toEqual([]); + }); + + it("caps the column count", () => { + expect(MAX_MODEL_FEATURES).toBe(12); + }); +}); + +describe("buildFeatureMatrix", () => { + function drivenDays(count: number): PrognosisDayInput[] { + return Array.from({ length: count }, (_, i) => + day(i, { a1: 10 - (i % 11), dimensions: { a2: i % 11 } }), + ); + } + + it("standardises per user, so the columns are that person's own units", () => { + const matrix = buildFeatureMatrix(drivenDays(40)); + expect(matrix).not.toBeNull(); + if (!matrix) return; + const column = matrix.rows.map((row) => row[0]); + const mean = column.reduce((s, v) => s + v, 0) / column.length; + const sd = Math.sqrt( + column.reduce((s, v) => s + (v - mean) * (v - mean), 0) / column.length, + ); + expect(mean).toBeCloseTo(0, 6); + expect(sd).toBeCloseTo(1, 6); + // The standardisation is kept so a later day can be scored on the same + // scale rather than re-standardised against itself. + expect(matrix.standardisation[0].key).toBe(matrix.features[0]); + expect(matrix.standardisation[0].sd).toBeGreaterThan(0); + }); + + it("orders rows oldest first, which is what makes the validation honest", () => { + const matrix = buildFeatureMatrix(drivenDays(40)); + expect(matrix).not.toBeNull(); + if (!matrix) return; + const sorted = [...matrix.days].sort(); + expect(matrix.days).toEqual(sorted); + }); + + it("drops a day that is missing a selected feature rather than filling it in", () => { + const days = drivenDays(40); + days[5] = { ...days[5], dimensions: {} }; + const matrix = buildFeatureMatrix(days); + expect(matrix).not.toBeNull(); + if (!matrix) return; + expect(matrix.days).not.toContain(days[5].day); + expect(matrix.rows).toHaveLength(39); + // Nothing was imputed in its place. + expect(matrix.targets).toHaveLength(39); + }); + + it("refuses when nothing survives the screening", () => { + const flat = Array.from({ length: 40 }, (_, i) => day(i, { a1: 5 })); + expect(buildFeatureMatrix(flat)).toBeNull(); + }); + + it("refuses a constant target instead of dividing by zero", () => { + const constant = Array.from({ length: 40 }, (_, i) => + day(i, { a1: 5, dimensions: { a2: i % 7 } }), + ); + expect(buildFeatureMatrix(constant)).toBeNull(); + }); +}); + +describe("standardiseDay", () => { + it("scores a day on the stored scale", () => { + const row = standardiseDay(day(0, { dimensions: { a2: 5 } }), [ + { key: "dimension:a2", mean: 5, sd: 2 }, + ]); + // Up-is-better A2 of a stress-5 day is 5, the mean, so zero. + expect(row).toEqual([0]); + }); + + it("refuses a day missing one of the fit's features", () => { + expect( + standardiseDay(day(0, { dimensions: {} }), [ + { key: "dimension:a2", mean: 5, sd: 2 }, + ]), + ).toBeNull(); + }); +}); diff --git a/src/lib/analytics/mood-prognosis/__tests__/rolling-origin.test.ts b/src/lib/analytics/mood-prognosis/__tests__/rolling-origin.test.ts new file mode 100644 index 000000000..53f93277a --- /dev/null +++ b/src/lib/analytics/mood-prognosis/__tests__/rolling-origin.test.ts @@ -0,0 +1,269 @@ +/** + * The fit, and the one property that makes its band worth reading. + * + * The leakage case is the reason this file exists. A random k-fold over days + * produces a model that scores beautifully and predicts nothing, and there is + * nothing in the reported error that says so — which is why the assertion is + * on the FOLD INDICES and on the day keys themselves rather than on the error. + * An error-based check passes a leaking split whenever the leak happens to be + * small, which is exactly when nobody would look. + * + * Watched red: the split was replaced with a shuffle, and + * `describe("no training fold contains a day after its validation fold")` + * went red naming a future day inside a training fold. Reverted, hash checked. + */ +import { describe, expect, it } from "vitest"; + +import { buildFeatureMatrix, type PrognosisDayInput } from "../features"; +import { fitRidge, predictWithFit, solveLinearSystem } from "../ridge"; +import { + BAND_LOWER_QUANTILE, + BAND_UPPER_QUANTILE, + LAMBDA_GRID, + MIN_VALIDATION_ROWS, + MODEL_VERSION, + fitPrognosisModel, + quantile, + rollingOriginFolds, +} from "../validate"; + +const EMPTY_LINKED = { + sleepAsleep: null, + steps: null, + activeEnergy: null, + restingHeartRate: null, + heartRateVariability: null, +}; + +function dayKey(index: number): string { + return new Date(Date.UTC(2026, 0, 1 + index)).toISOString().slice(0, 10); +} + +/** A deterministic pseudo-random series — a flaky statistics test is useless. */ +function seeded(seed: number): () => number { + let state = seed; + return () => { + state = (state * 1103515245 + 12345) % 2147483648; + return state / 2147483648; + }; +} + +/** + * Days whose target is a clean linear function of one dimension. + * + * One driver rather than two, and that is a statement about the screening + * rather than about the fit: the screening is univariate, so a feature that + * only matters jointly with another is refused before the fit ever sees it. + * A two-driver fixture measures that property, not this one. + */ +function linearDays(count: number): PrognosisDayInput[] { + return Array.from({ length: count }, (_, i) => { + const stress = i % 11; + return { + day: dayKey(i), + // 0.8 × up-is-better stress, plus one — inside the 0-10 scale. + a1: 0.8 * (10 - stress) + 1, + dimensions: { a2: stress, a3: (i * 7) % 11 }, + context: null, + linked: { ...EMPTY_LINKED }, + }; + }); +} + +/** Days whose target is noise and whose inputs are unrelated noise. */ +function noisyDays(count: number, seed = 7): PrognosisDayInput[] { + const random = seeded(seed); + return Array.from({ length: count }, (_, i) => ({ + day: dayKey(i), + a1: Math.round(random() * 10), + dimensions: { + a2: Math.round(random() * 10), + a3: Math.round(random() * 10), + }, + context: null, + linked: { ...EMPTY_LINKED }, + })); +} + +describe("no training fold contains a day after its validation fold", () => { + const ROWS = 60; + const folds = rollingOriginFolds(ROWS); + + it("produces folds at all — the assertions below pass on an empty list", () => { + expect(folds.length).toBeGreaterThan(0); + expect(folds.length).toBeLessThanOrEqual(4); + // Every validation row is covered exactly once and the last fold reaches + // the end, so no day is silently dropped from the measurement. + expect(folds[folds.length - 1].validationEnd).toBe(ROWS); + }); + + it("trains on a prefix that ends exactly where the validation starts", () => { + for (const fold of folds) { + expect(fold.trainEnd).toBe(fold.validationStart); + expect(fold.trainEnd).toBeGreaterThan(0); + expect(fold.validationEnd).toBeGreaterThan(fold.validationStart); + } + }); + + it("holds on the day keys, not just on the indices", () => { + // The indices are only meaningful because the matrix is time-ordered, so + // the property is asserted against the dates the rows actually carry. + const matrix = buildFeatureMatrix(linearDays(ROWS)); + expect(matrix).not.toBeNull(); + if (!matrix) return; + + for (const fold of rollingOriginFolds(matrix.rows.length)) { + const trainingDays = matrix.days.slice(0, fold.trainEnd); + const validationDays = matrix.days.slice( + fold.validationStart, + fold.validationEnd, + ); + const latestTrainingDay = trainingDays[trainingDays.length - 1]; + const earliestValidationDay = validationDays[0]; + const leaked = trainingDays.filter((d) => d >= earliestValidationDay); + expect( + leaked, + `a training fold holds ${leaked.length} day(s) dated at or after the validation fold that starts on ${earliestValidationDay}`, + ).toEqual([]); + expect(latestTrainingDay < earliestValidationDay).toBe(true); + } + }); + + it("expands rather than slides — later folds train on more history", () => { + const sizes = folds.map((f) => f.trainEnd); + for (let i = 1; i < sizes.length; i++) { + expect(sizes[i]).toBeGreaterThan(sizes[i - 1]); + } + }); + + it("refuses to split a window too short to measure anything on", () => { + expect(rollingOriginFolds(MIN_VALIDATION_ROWS - 1)).toEqual([]); + }); +}); + +describe("the ridge solve", () => { + it("recovers the coefficients of a clean linear signal", () => { + // y = 2*x1 - 3*x2, standardised columns, tiny penalty. + const rows: number[][] = []; + const targets: number[] = []; + for (let i = 0; i < 40; i++) { + const x1 = ((i % 7) - 3) / 2; + const x2 = ((i % 5) - 2) / 1.5; + rows.push([x1, x2]); + targets.push(2 * x1 - 3 * x2 + 5); + } + const fit = fitRidge(rows, targets, 1e-6); + expect(fit).not.toBeNull(); + if (!fit) return; + expect(fit.coefficients[0]).toBeCloseTo(2, 2); + expect(fit.coefficients[1]).toBeCloseTo(-3, 2); + // The intercept is the mean of the training targets, and the columns are + // centred against their own means, so a row at the column means predicts + // the average day exactly. + expect(predictWithFit(fit, fit.columnMeans)).toBeCloseTo(fit.intercept, 6); + // And every training row is predicted back to its own target. + rows.forEach((row, i) => { + expect(predictWithFit(fit, row)).toBeCloseTo(targets[i], 2); + }); + }); + + it("does not divide by zero on a constant target", () => { + const rows = Array.from({ length: 20 }, (_, i) => [i / 10, (i % 3) / 2]); + const fit = fitRidge(rows, new Array(20).fill(4), 1); + expect(fit).not.toBeNull(); + if (!fit) return; + expect(fit.coefficients.every(Number.isFinite)).toBe(true); + expect(fit.intercept).toBe(4); + expect(predictWithFit(fit, [0, 0])).toBeCloseTo(4, 6); + }); + + it("refuses a singular system rather than answering zeros", () => { + // Two identical columns and no penalty: the system has no unique answer, + // and a vector of zeros would look exactly like "these features do not + // matter". + expect( + solveLinearSystem( + [ + [1, 1], + [1, 1], + ], + [1, 1], + ), + ).toBeNull(); + }); + + it("penalises: a larger lambda pulls the coefficients towards zero", () => { + const rows = Array.from({ length: 30 }, (_, i) => [((i % 7) - 3) / 2]); + const targets = rows.map((r) => 4 * r[0] + 5); + const light = fitRidge(rows, targets, 0.01); + const heavy = fitRidge(rows, targets, 100); + expect(light).not.toBeNull(); + expect(heavy).not.toBeNull(); + if (!light || !heavy) return; + expect(Math.abs(heavy.coefficients[0])).toBeLessThan( + Math.abs(light.coefficients[0]), + ); + }); +}); + +describe("the band reacts to how wrong the model is", () => { + it("is narrow on a clean signal and wide on noise", () => { + const clean = buildFeatureMatrix(linearDays(80)); + expect(clean).not.toBeNull(); + if (!clean) return; + const cleanModel = fitPrognosisModel(clean); + expect(cleanModel).not.toBeNull(); + if (!cleanModel) return; + const cleanWidth = cleanModel.bandHigh - cleanModel.bandLow; + + // Noise cannot go through the screening — it is refused before a fit is + // attempted, which is the stronger answer. So the band's reaction is + // measured on the matrix directly, with the screening bypassed by handing + // the fit a matrix built from the noisy targets and the clean columns. + const noisy = buildFeatureMatrix(linearDays(80)); + expect(noisy).not.toBeNull(); + if (!noisy) return; + const random = seeded(11); + const noisyMatrix = { + ...noisy, + targets: noisy.targets.map(() => Math.round(random() * 10)), + }; + const noisyModel = fitPrognosisModel(noisyMatrix); + expect(noisyModel).not.toBeNull(); + if (!noisyModel) return; + const noisyWidth = noisyModel.bandHigh - noisyModel.bandLow; + + expect( + noisyWidth, + "the band did not widen on noise — a confident wrong number is the failure this band exists to prevent", + ).toBeGreaterThan(cleanWidth * 2); + expect(cleanWidth).toBeLessThan(2); + }); + + it("refuses noise at the screening, before any fit is attempted", () => { + expect(buildFeatureMatrix(noisyDays(80))).toBeNull(); + }); + + it("takes the band from the validation residuals, at the stated quantiles", () => { + const matrix = buildFeatureMatrix(linearDays(80)); + expect(matrix).not.toBeNull(); + if (!matrix) return; + const model = fitPrognosisModel(matrix); + expect(model).not.toBeNull(); + if (!model) return; + expect(model.bandLow).toBeLessThanOrEqual(model.bandHigh); + expect(BAND_LOWER_QUANTILE).toBe(0.1); + expect(BAND_UPPER_QUANTILE).toBe(0.9); + expect(model.n).toBe(matrix.rows.length); + expect(model.modelVersion).toBe(MODEL_VERSION); + expect(LAMBDA_GRID).toContain(model.fit.lambda); + // Never zero: an unpenalised solve memorises the training days. + expect(Math.min(...LAMBDA_GRID)).toBeGreaterThan(0); + }); + + it("quantile interpolates rather than rounding to a sample", () => { + expect(quantile([0, 10], 0.5)).toBe(5); + expect(quantile([1, 2, 3, 4], 0)).toBe(1); + expect(quantile([1, 2, 3, 4], 1)).toBe(4); + }); +}); diff --git a/src/lib/analytics/mood-prognosis/feature-keys.ts b/src/lib/analytics/mood-prognosis/feature-keys.ts new file mode 100644 index 000000000..3f51c71b5 --- /dev/null +++ b/src/lib/analytics/mood-prognosis/feature-keys.ts @@ -0,0 +1,98 @@ +/** + * The names the model's inputs are known by, everywhere. + * + * One string per feature, stored in `MoodPrediction.features`, read back by + * the surface that names what the model weighted, and compared against by the + * tests. It is therefore a wire format in all but name: renaming a key + * silently orphans every stored row that used it, which is why the encoding is + * a function here rather than a template literal at four call sites. + * + * The shape is `:`, and the four families are the four places an + * input can come from: + * + * * `dimension:a2` — a level-A self-rating (A1 is the target) + * * `context:workLoad` — a numeric field of the day's context + * * `context:workStatus=off` — one value of a closed context vocabulary + * * `linked:steps` — a figure read from the module that owns it + * + * Nothing here resolves a label. The label is a locale lookup and this module + * is pure and bundle-safe, so the resolver lives beside the component that + * renders it and takes a parsed key from here. + */ +import type { MoodDimensionKey } from "@/lib/mood/dimensions"; + +/** The linked figures the model may read, and the day-figure each names. */ +export const LINKED_FEATURE_METRICS = [ + "sleepAsleep", + "steps", + "activeEnergy", + "restingHeartRate", + "heartRateVariability", +] as const; + +export type LinkedFeatureMetric = (typeof LINKED_FEATURE_METRICS)[number]; + +export type ParsedFeatureKey = + | { family: "dimension"; dimension: MoodDimensionKey } + | { family: "context"; field: string; value: null } + | { family: "context"; field: string; value: string } + | { family: "linked"; metric: LinkedFeatureMetric }; + +/** `dimension:a2` */ +export function dimensionFeatureKey(dimension: MoodDimensionKey): string { + return `dimension:${dimension}`; +} + +/** `context:workLoad` — a numeric field. */ +export function contextNumericFeatureKey(field: string): string { + return `context:${field}`; +} + +/** `context:workStatus=off` — one value of a closed vocabulary. */ +export function contextValueFeatureKey(field: string, value: string): string { + return `context:${field}=${value}`; +} + +/** `linked:steps` */ +export function linkedFeatureKey(metric: LinkedFeatureMetric): string { + return `linked:${metric}`; +} + +const DIMENSION_KEYS = new Set(["a1", "a2", "a3", "a4", "a5"]); +const LINKED_KEYS = new Set(LINKED_FEATURE_METRICS); + +/** + * Read a stored key back. + * + * `null` for anything this version does not recognise, and that case is real + * rather than defensive: a row written by an older `modelVersion` can name a + * feature the current build has dropped, and the surface has to skip it rather + * than render `undefined` at somebody. + */ +export function parseFeatureKey(key: string): ParsedFeatureKey | null { + const separator = key.indexOf(":"); + if (separator < 0) return null; + const family = key.slice(0, separator); + const rest = key.slice(separator + 1); + if (rest.length === 0) return null; + + if (family === "dimension") { + return DIMENSION_KEYS.has(rest) + ? { family: "dimension", dimension: rest as MoodDimensionKey } + : null; + } + if (family === "linked") { + return LINKED_KEYS.has(rest) + ? { family: "linked", metric: rest as LinkedFeatureMetric } + : null; + } + if (family === "context") { + const equals = rest.indexOf("="); + if (equals < 0) return { family: "context", field: rest, value: null }; + const field = rest.slice(0, equals); + const value = rest.slice(equals + 1); + if (field.length === 0 || value.length === 0) return null; + return { family: "context", field, value }; + } + return null; +} diff --git a/src/lib/analytics/mood-prognosis/features.ts b/src/lib/analytics/mood-prognosis/features.ts new file mode 100644 index 000000000..467b71236 --- /dev/null +++ b/src/lib/analytics/mood-prognosis/features.ts @@ -0,0 +1,489 @@ +/** + * One row per day, one column per thing the day was made of. + * + * The target is the day's own pleasantness (A1) and everything else is an + * input: the other four level-A dimensions, the numbers and the closed + * vocabularies of the day's context, and the figures the sleep, activity and + * recovery modules already hold for that day. A1 is the only target on + * purpose — five targets would be five ways to be wrong and five bands to + * explain. + * + * Four rules this file is built on, each of which has a way of being got + * wrong that looks fine on screen: + * + * 1. **Nothing is imputed. Ever.** A day with a missing input is not a day + * with an average input. The stated rule is complete cases: a day enters + * the fit only if it carries the target and every SELECTED feature, and a + * feature only becomes a candidate if it was recorded on at least + * {@link MIN_FEATURE_COVERAGE} of the days that carry a target. The + * coverage floor is what keeps the complete-case rule from quietly + * emptying the matrix: without it, one field somebody filled in twice + * would exclude every other day they ever logged. + * 2. **Standardised per user, not min-max.** A person whose stress never + * exceeds six must not have six read as "maximum stress". The mean and the + * spread come from their own fit window, so every feature is expressed in + * that person's own units — which is also what makes the stored + * contributions comparable to each other. + * 3. **Screened before it is fitted, with the multiple-comparison correction + * the rest of the product already runs.** A sweep over sixty context + * features at raw p < 0.05 finds three things that are not there. Pearson + * with the shared `MIN_PAIRED_N` floor, Benjamini-Hochberg at the shared + * `FDR_Q`, and the shared effect-size floor on the shrunk estimate — all + * imported, none re-tuned here. + * 4. **Orientation comes from the table that owns it.** A2 is inverse and is + * flipped through `MOOD_DIMENSIONS[i].inverse`, never inline, so the chart + * and the model read one source. Context ratings are NOT flipped: no table + * states their orientation, and inventing one here would be exactly the + * second opinion the dimension table exists to prevent. It costs nothing + * mathematically — a flip only changes the sign of a coefficient, and the + * stored contribution is signed anyway. + */ +import { + EFFECT_SIZE_FLOOR, + FDR_Q, + benjaminiHochberg, + shrinkEstimate, +} from "@/lib/insights/correlation-discovery"; +import { + MAX_P_VALUE, + MIN_PAIRED_N, + pearson, +} from "@/lib/insights/correlations"; +import { + CONTACT_CIRCLE_KEYS, + CONTACT_EXTENT_KEYS, + CONTACT_FORM_KEYS, + EVENT_TYPE_KEYS, + LEISURE_CATEGORY_KEYS, + WORK_STATUS_KEYS, +} from "@/lib/mood/context-vocabulary"; +import { MOOD_DIMENSIONS } from "@/lib/mood/dimensions"; +import type { MoodDimensionKey } from "@/lib/mood/dimensions"; + +import { + contextNumericFeatureKey, + contextValueFeatureKey, + dimensionFeatureKey, + linkedFeatureKey, + LINKED_FEATURE_METRICS, +} from "./feature-keys"; +import { occurrencesSuffice } from "./thresholds"; + +/** + * Most features the fit may carry. + * + * Twelve over thirty days is already generous — a ridge penalty is what makes + * that survivable at all — and the cap is what stops a well-logged account + * from being fitted with sixty columns and a hundred rows. + */ +export const MAX_MODEL_FEATURES = 12; + +/** + * Share of target-carrying days a feature must be recorded on to be a + * candidate at all. + * + * The other half of the complete-case rule. A field somebody filled in for a + * fortnight and then stopped is not a column of the matrix; it is a fortnight, + * and admitting it would silently discard every day outside that fortnight. + */ +export const MIN_FEATURE_COVERAGE = 0.6; + +/** The day's context as the model reads it. The note never comes here. */ +export interface PrognosisContextInput { + workStatus: string | null; + workMinutes: number | null; + overtimeMinutes: number | null; + workLoad: number | null; + workSatisfaction: number | null; + contactCircles: string[]; + contactForm: string | null; + contactExtent: string | null; + contactQuality: number | null; + contactSupport: number | null; + leisureCategories: string[]; + leisureMinutes: number | null; + leisureJoy: number | null; + leisureRecovery: number | null; + eventType: string | null; + eventValence: number | null; +} + +/** One day, as everything upstream of the fit has resolved it. */ +export interface PrognosisDayInput { + /** Local day key, `YYYY-MM-DD`. */ + day: string; + /** The target: the day's own pleasantness, 0-10. `null` days are not fitted. */ + a1: number | null; + /** The other four self-ratings, literal as the person set them. */ + dimensions: Partial>; + /** The day's stored context, or `null` when no section was opened. */ + context: PrognosisContextInput | null; + /** What the owning modules hold for the day. `null` is absence. */ + linked: { + sleepAsleep: number | null; + steps: number | null; + activeEnergy: number | null; + restingHeartRate: number | null; + heartRateVariability: number | null; + }; +} + +/** The numeric context fields, in capture order. */ +const CONTEXT_NUMERIC_FIELDS = [ + "workMinutes", + "overtimeMinutes", + "workLoad", + "workSatisfaction", + "contactQuality", + "contactSupport", + "leisureMinutes", + "leisureJoy", + "leisureRecovery", + "eventValence", +] as const satisfies readonly (keyof PrognosisContextInput)[]; + +/** The single-valued vocabularies, each expanded to one column per value. */ +const CONTEXT_SINGLE_VOCABULARIES = [ + { field: "workStatus", keys: WORK_STATUS_KEYS }, + { field: "contactForm", keys: CONTACT_FORM_KEYS }, + { field: "contactExtent", keys: CONTACT_EXTENT_KEYS }, + { field: "eventType", keys: EVENT_TYPE_KEYS }, +] as const; + +/** The multi-select vocabularies, same expansion. */ +const CONTEXT_MULTI_VOCABULARIES = [ + { field: "contactCircles", keys: CONTACT_CIRCLE_KEYS }, + { field: "leisureCategories", keys: LEISURE_CATEGORY_KEYS }, +] as const; + +function isNumber(value: unknown): value is number { + return typeof value === "number" && Number.isFinite(value); +} + +/** + * Every feature this day carries a value for. + * + * A feature the day says nothing about is ABSENT from the map — it is not + * zero. The difference is the whole of rule 1: a day with no context row has + * no work-status column at all, while a day whose work status is `off` has a + * zero in every other work-status column and a one in that one. + */ +export function extractDayFeatures( + day: PrognosisDayInput, +): Map { + const out = new Map(); + + for (const dimension of MOOD_DIMENSIONS) { + // A1 is the target, not an input. Fitting a value on itself would produce + // a perfect model of nothing. + if (dimension.key === "a1") continue; + const raw = day.dimensions[dimension.key]; + if (!isNumber(raw)) continue; + // The orientation flip goes through the table that declares it, so a + // dimension that becomes inverse later moves the chart and the model + // together. + out.set( + dimensionFeatureKey(dimension.key), + dimension.inverse ? dimension.max - raw : raw, + ); + } + + const context = day.context; + if (context) { + for (const field of CONTEXT_NUMERIC_FIELDS) { + const raw = context[field]; + if (isNumber(raw)) out.set(contextNumericFeatureKey(field), raw); + } + for (const { field, keys } of CONTEXT_SINGLE_VOCABULARIES) { + const value = context[field]; + // A field nobody answered contributes no columns. Zeroing them would + // tell the fit "not overtime" on a day nobody said anything about work. + if (typeof value !== "string" || value.length === 0) continue; + for (const key of keys) { + out.set(contextValueFeatureKey(field, key), value === key ? 1 : 0); + } + } + for (const { field, keys } of CONTEXT_MULTI_VOCABULARIES) { + const values = context[field]; + if (!Array.isArray(values) || values.length === 0) continue; + const chosen = new Set(values); + for (const key of keys) { + out.set(contextValueFeatureKey(field, key), chosen.has(key) ? 1 : 0); + } + } + } + + for (const metric of LINKED_FEATURE_METRICS) { + const raw = day.linked[metric]; + if (isNumber(raw)) out.set(linkedFeatureKey(metric), raw); + } + + return out; +} + +/** What the screening decided about one candidate, kept for the tests. */ +export interface ScreenedFeature { + key: string; + /** Days carrying both the feature and the target. */ + pairs: number; + r: number; + /** The shrunk estimate the effect floor is applied to. */ + shrunk: number; + pValue: number; + qValue: number; + selected: boolean; + /** Why it was refused, when it was. */ + reason?: + | "coverage" + | "too-few-pairs" + | "no-variance" + | "p-value" + | "fdr" + | "effect-size" + | "occurrences" + | "cap"; +} + +/** + * Which candidates earn a column, and why the others did not. + * + * Order of refusal matters and is the order below: coverage, then pairs, then + * the raw p, then the family-wide FDR, then the effect floor, then the cap. + * A feature that fails an earlier gate is not tested against a later one, so + * the family the FDR corrects over is the set that was actually tested. + */ +export function screenFeatures( + days: readonly PrognosisDayInput[], +): ScreenedFeature[] { + const targetDays = days.filter((d) => isNumber(d.a1)); + const perDay = new Map>(); + for (const day of targetDays) perDay.set(day.day, extractDayFeatures(day)); + + const candidates = new Set(); + for (const features of perDay.values()) { + for (const key of features.keys()) candidates.add(key); + } + + const coverageFloor = targetDays.length * MIN_FEATURE_COVERAGE; + const refused: ScreenedFeature[] = []; + const tested: Array<{ + key: string; + pairs: number; + r: number; + shrunk: number; + pValue: number; + }> = []; + + for (const key of [...candidates].sort()) { + const xs: number[] = []; + const ys: number[] = []; + let occurrences = 0; + for (const day of targetDays) { + const value = perDay.get(day.day)?.get(key); + if (value === undefined) continue; + xs.push(value); + ys.push(day.a1 as number); + if (value !== 0) occurrences += 1; + } + + if (xs.length < coverageFloor) { + refused.push(zero(key, xs.length, "coverage")); + continue; + } + // A one-hot column seen four times is not evidence, however well it + // correlates. This is the ladder's occurrence floor, not a second one. + if (!occurrencesSuffice(occurrences)) { + refused.push(zero(key, xs.length, "occurrences")); + continue; + } + + const result = pearson({ xs, ys, minPairs: MIN_PAIRED_N }); + if (result.status === "insufficient") { + refused.push( + zero( + key, + result.n, + result.reason === "no_variance" ? "no-variance" : "too-few-pairs", + ), + ); + continue; + } + if (result.pValue >= MAX_P_VALUE) { + refused.push({ + key, + pairs: result.n, + r: result.r, + shrunk: shrinkEstimate(result.r, result.n), + pValue: result.pValue, + qValue: 1, + selected: false, + reason: "p-value", + }); + continue; + } + tested.push({ + key, + pairs: result.n, + r: result.r, + shrunk: shrinkEstimate(result.r, result.n), + pValue: result.pValue, + }); + } + + const qValues = benjaminiHochberg(tested.map((t) => t.pValue)); + const survivors: ScreenedFeature[] = []; + tested.forEach((t, index) => { + const qValue = qValues[index] ?? 1; + if (qValue > FDR_Q) { + survivors.push({ ...t, qValue, selected: false, reason: "fdr" }); + return; + } + if (Math.abs(t.shrunk) < EFFECT_SIZE_FLOOR) { + survivors.push({ ...t, qValue, selected: false, reason: "effect-size" }); + return; + } + survivors.push({ ...t, qValue, selected: true }); + }); + + // The cap keeps the strongest, and a stable tie-break by key keeps two runs + // over the same data returning the same columns. + const admitted = survivors + .filter((s) => s.selected) + .sort( + (a, b) => + Math.abs(b.shrunk) - Math.abs(a.shrunk) || a.key.localeCompare(b.key), + ); + for (const feature of admitted.slice(MAX_MODEL_FEATURES)) { + feature.selected = false; + feature.reason = "cap"; + } + + return [...survivors, ...refused].sort((a, b) => a.key.localeCompare(b.key)); +} + +function zero( + key: string, + pairs: number, + reason: ScreenedFeature["reason"], +): ScreenedFeature { + return { + key, + pairs, + r: 0, + shrunk: 0, + pValue: 1, + qValue: 1, + selected: false, + reason, + }; +} + +/** Per-feature centre and spread, over this account's own fit window. */ +export interface FeatureStandardisation { + key: string; + mean: number; + sd: number; +} + +export interface FeatureMatrix { + /** The selected feature keys, in column order. */ + features: string[]; + /** Standardised rows, one per complete day, oldest first. */ + rows: number[][]; + /** The target for each row, unstandardised (0-10). */ + targets: number[]; + /** The day each row is, so the validation can stay time-ordered. */ + days: string[]; + /** What each column was standardised by, needed to score a new day. */ + standardisation: FeatureStandardisation[]; +} + +/** + * Build the matrix: screen, then keep the complete days, then standardise. + * + * Returns `null` when nothing survives the screening or when no day carries + * every surviving feature. That is a refusal, not an empty matrix — the caller + * writes no forecast rather than one built on nothing. + */ +export function buildFeatureMatrix( + days: readonly PrognosisDayInput[], +): FeatureMatrix | null { + const screened = screenFeatures(days); + const selected = screened.filter((s) => s.selected).map((s) => s.key); + if (selected.length === 0) return null; + + const ordered = [...days] + .filter((d) => isNumber(d.a1)) + .sort((a, b) => a.day.localeCompare(b.day)); + + const rawRows: number[][] = []; + const targets: number[] = []; + const dayKeys: string[] = []; + for (const day of ordered) { + const features = extractDayFeatures(day); + const row: number[] = []; + let complete = true; + for (const key of selected) { + const value = features.get(key); + if (value === undefined) { + complete = false; + break; + } + row.push(value); + } + // Complete cases only. A day missing one selected feature is a day the fit + // does not see; nothing is filled in for it. + if (!complete) continue; + rawRows.push(row); + targets.push(day.a1 as number); + dayKeys.push(day.day); + } + + if (rawRows.length === 0) return null; + + const standardisation = selected.map((key, column) => { + const values = rawRows.map((row) => row[column]); + const mean = values.reduce((sum, v) => sum + v, 0) / values.length; + const variance = + values.reduce((sum, v) => sum + (v - mean) * (v - mean), 0) / + values.length; + return { key, mean, sd: Math.sqrt(variance) }; + }); + + // A column with no spread inside the complete-case set carries no + // information and would divide by zero. It goes, and if that empties the + // set the whole fit is refused. + const keep = standardisation + .map((s, index) => ({ ...s, index })) + .filter((s) => s.sd > 0); + if (keep.length === 0) return null; + + return { + features: keep.map((k) => k.key), + rows: rawRows.map((row) => keep.map((k) => (row[k.index] - k.mean) / k.sd)), + targets, + days: dayKeys, + standardisation: keep.map(({ key, mean, sd }) => ({ key, mean, sd })), + }; +} + +/** + * Standardise one day's raw features against a stored fit, for scoring. + * + * `null` when the day does not carry every feature the fit uses — the same + * complete-case rule the matrix applies, asked of a single day. A day that + * cannot be scored gets no forecast, which is the honest answer. + */ +export function standardiseDay( + day: PrognosisDayInput, + standardisation: readonly FeatureStandardisation[], +): number[] | null { + const features = extractDayFeatures(day); + const row: number[] = []; + for (const { key, mean, sd } of standardisation) { + const value = features.get(key); + if (value === undefined || sd <= 0) return null; + row.push((value - mean) / sd); + } + return row; +} diff --git a/src/lib/analytics/mood-prognosis/ridge.ts b/src/lib/analytics/mood-prognosis/ridge.ts new file mode 100644 index 000000000..8f5335a4f --- /dev/null +++ b/src/lib/analytics/mood-prognosis/ridge.ts @@ -0,0 +1,173 @@ +/** + * Ridge regression, closed form, in about sixty lines of arithmetic. + * + * `(XᵀX + λI)β = Xᵀy`, solved by Gaussian elimination with partial pivoting. + * That is the whole model. It is here rather than in a dependency because a + * dozen standardised columns is a 12×12 system, and adding a package to solve + * one of those would be adding a supply chain to avoid writing a loop. + * + * Two properties the rest of the feature relies on: + * + * * **The intercept is not penalised, because there is no intercept in the + * system.** The target is centred before the solve and the mean is added + * back afterwards, which is the standard treatment: penalising the + * intercept would shrink every prediction towards zero rather than + * towards the person's own average day, and towards zero on a 0-10 + * pleasantness scale means "towards the worst day imaginable". + * * **λ is never zero.** With more columns than the ridge penalty can + * afford, an unpenalised solve is a memoriser: it fits the training days + * perfectly and says nothing about any other day. The grid the validation + * searches starts above zero for that reason. + * + * This is arithmetic, not an AI surface. No provider, no consent gate, no + * egress — see the note at the top of the refresh job. + */ + +/** A fitted model: one coefficient per standardised column, plus the mean. */ +export interface RidgeFit { + /** One per column of the matrix it was fitted on, in the same order. */ + coefficients: number[]; + /** The mean of the training targets, added back on every prediction. */ + intercept: number; + /** + * The mean of each training column, subtracted before every prediction. + * + * Usually zero to floating point, because the matrix arrives z-scored. It is + * carried anyway rather than assumed: centring the target while leaving the + * columns uncentred biases every coefficient, and the bias is small enough + * to look like rounding. Somebody handing this function a raw matrix should + * get the right answer, not a plausible one. + */ + columnMeans: number[]; + /** The penalty this fit used. */ + lambda: number; +} + +/** + * Solve `A x = b` by Gaussian elimination with partial pivoting. + * + * `null` when the system is singular to working precision. The caller treats + * that as "no fit" rather than as a fit of zeros: a singular system means the + * columns are telling the same story twice, and a zero vector would silently + * become a model that predicts the average and calls it a pattern. + */ +export function solveLinearSystem( + matrix: readonly number[][], + vector: readonly number[], +): number[] | null { + const n = vector.length; + if (matrix.length !== n) return null; + // Copy: the caller's XᵀX is reused across the λ grid and elimination is + // destructive. + const a = matrix.map((row, i) => [...row, vector[i]]); + + for (let column = 0; column < n; column++) { + let pivot = column; + for (let row = column + 1; row < n; row++) { + if (Math.abs(a[row][column]) > Math.abs(a[pivot][column])) pivot = row; + } + if (Math.abs(a[pivot][column]) < 1e-10) return null; + if (pivot !== column) { + const swap = a[pivot]; + a[pivot] = a[column]; + a[column] = swap; + } + const pivotValue = a[column][column]; + for (let row = column + 1; row < n; row++) { + const factor = a[row][column] / pivotValue; + if (factor === 0) continue; + for (let k = column; k <= n; k++) a[row][k] -= factor * a[column][k]; + } + } + + const solution = new Array(n).fill(0); + for (let row = n - 1; row >= 0; row--) { + let sum = a[row][n]; + for (let k = row + 1; k < n; k++) sum -= a[row][k] * solution[k]; + solution[row] = sum / a[row][row]; + } + return solution.every((v) => Number.isFinite(v)) ? solution : null; +} + +/** + * Fit `rows` against `targets` at one penalty. + * + * `rows` are already standardised per user (see `features.ts`), so the columns + * are comparable and one λ is a defensible penalty across all of them — + * which is not true of raw minutes beside a 0-10 rating. + */ +export function fitRidge( + rows: readonly (readonly number[])[], + targets: readonly number[], + lambda: number, +): RidgeFit | null { + const n = rows.length; + if (n === 0 || targets.length !== n) return null; + const columns = rows[0].length; + if (columns === 0) return null; + if (rows.some((row) => row.length !== columns)) return null; + + const intercept = targets.reduce((sum, v) => sum + v, 0) / n; + const centredTargets = targets.map((v) => v - intercept); + const columnMeans = Array.from({ length: columns }, (_, j) => { + let sum = 0; + for (let i = 0; i < n; i++) sum += rows[i][j]; + return sum / n; + }); + const centredRows = rows.map((row) => row.map((v, j) => v - columnMeans[j])); + + const xtx: number[][] = Array.from({ length: columns }, () => + new Array(columns).fill(0), + ); + const xty = new Array(columns).fill(0); + for (let i = 0; i < n; i++) { + const row = centredRows[i]; + for (let j = 0; j < columns; j++) { + xty[j] += row[j] * centredTargets[i]; + for (let k = j; k < columns; k++) { + const product = row[j] * row[k]; + xtx[j][k] += product; + if (k !== j) xtx[k][j] += product; + } + } + } + for (let j = 0; j < columns; j++) xtx[j][j] += lambda; + + const coefficients = solveLinearSystem(xtx, xty); + if (coefficients === null) return null; + return { coefficients, intercept, columnMeans, lambda }; +} + +/** Score one standardised row. */ +export function predictWithFit( + fit: RidgeFit, + row: readonly number[], +): number | null { + if (row.length !== fit.coefficients.length) return null; + let sum = fit.intercept; + for (let i = 0; i < row.length; i++) { + sum += fit.coefficients[i] * (row[i] - fit.columnMeans[i]); + } + return Number.isFinite(sum) ? sum : null; +} + +/** + * The per-feature contributions of one scored day. + * + * `coefficient × the day's standardised value`, which is what the stored + * `features` column holds and what the explanation line is generated from. + * Signed, so a feature that pulled the expected value down reads as having + * pulled it down, and the parts sum to the distance between the prediction and + * the person's own average day. + */ +export function featureContributions( + fit: RidgeFit, + featureKeys: readonly string[], + row: readonly number[], +): Array<{ feature: string; contribution: number }> { + return featureKeys.map((feature, index) => ({ + feature, + contribution: + fit.coefficients[index] * (row[index] - fit.columnMeans[index]), + })); +} diff --git a/src/lib/analytics/mood-prognosis/validate.ts b/src/lib/analytics/mood-prognosis/validate.ts new file mode 100644 index 000000000..d5568c2e4 --- /dev/null +++ b/src/lib/analytics/mood-prognosis/validate.ts @@ -0,0 +1,202 @@ +/** + * Choosing the penalty, and measuring how wrong the model is, without letting + * it see the future. + * + * **Rolling origin, expanding window. Never a random k-fold.** Days are not + * exchangeable: a Tuesday sits between a Monday and a Wednesday that resemble + * it, and a shuffled split puts Wednesday in the training set and asks the + * model about Tuesday. The result is a model that looks excellent and is + * useless, and nothing about the number it reports says so. So the folds are + * ranges of a time-ordered array, the training range always ends where the + * validation range starts, and `rolling-origin.test.ts` asserts that on the + * fold indices and on the day keys themselves rather than on the resulting + * error — an error-based assertion would pass a leaking split whenever the + * leak happened to be small. + * + * **The band comes from the validation residuals**, not from a textbook + * standard error on the training fit. A training-fit interval on a model with + * twelve columns and thirty rows is narrow by construction and says nothing; + * the spread of the errors the model actually made on days it had not seen is + * the honest width. It widens when the model is bad, which is the property the + * surface depends on to report low confidence instead of a confident wrong + * number. + */ +import { fitRidge, predictWithFit, type RidgeFit } from "./ridge"; +import type { FeatureMatrix } from "./features"; + +/** + * The stored model identity. Bump it whenever the feature set, the screening + * or the fit changes, so a re-fit is legible in the rows rather than silently + * changing what they mean. Rows carrying an older version stay readable; the + * surface skips a feature key the current build no longer knows. + */ +export const MODEL_VERSION = "mood-ridge-1"; + +/** + * The penalties the validation searches. + * + * Never zero: with a dozen standardised columns over a few dozen days, an + * unpenalised solve memorises the training days and describes no other day. + * The grid is coarse on purpose — the difference between λ=8 and λ=10 is far + * below the noise in thirty days of self-ratings, and a fine grid would only + * be a slower way to overfit the validation set. + */ +export const LAMBDA_GRID = [0.5, 1, 2, 5, 10, 25, 50] as const; + +/** How many validation blocks the window is cut into. */ +export const VALIDATION_FOLDS = 4; + +/** + * The band is the central 80% of the validation residuals: the 10th and the + * 90th percentile. A range wide enough to be honest and narrow enough to be + * worth reading — and empirical, so an asymmetric error stays asymmetric + * instead of being averaged into a symmetric claim. + */ +export const BAND_LOWER_QUANTILE = 0.1; +export const BAND_UPPER_QUANTILE = 0.9; + +/** Fewest rows a fold split is attempted on at all. */ +export const MIN_VALIDATION_ROWS = 10; + +/** One expanding-window split. Every index is into the time-ordered matrix. */ +export interface RollingOriginFold { + /** Training rows are `[0, trainEnd)`. */ + trainEnd: number; + /** Validation rows are `[validationStart, validationEnd)`. */ + validationStart: number; + validationEnd: number; +} + +/** + * Cut a time-ordered window into expanding-window folds. + * + * `trainEnd === validationStart` in every fold, which is the whole invariant: + * the training range is a prefix and the validation range starts exactly where + * it ends. No gap (a gap would waste days), no overlap (an overlap is the + * leak). + */ +export function rollingOriginFolds( + rowCount: number, + folds: number = VALIDATION_FOLDS, +): RollingOriginFold[] { + if (rowCount < MIN_VALIDATION_ROWS || folds < 1) return []; + const block = Math.floor(rowCount / (folds + 1)); + if (block < 1) return []; + + const out: RollingOriginFold[] = []; + for (let index = 1; index <= folds; index++) { + const validationStart = block * index; + const validationEnd = + index === folds ? rowCount : Math.min(rowCount, block * (index + 1)); + if (validationEnd <= validationStart) continue; + out.push({ trainEnd: validationStart, validationStart, validationEnd }); + } + return out; +} + +/** Linear-interpolated quantile of an unsorted sample. */ +export function quantile(values: readonly number[], q: number): number { + if (values.length === 0) return 0; + const sorted = [...values].sort((a, b) => a - b); + if (sorted.length === 1) return sorted[0]; + const position = (sorted.length - 1) * Math.min(1, Math.max(0, q)); + const lower = Math.floor(position); + const upper = Math.ceil(position); + if (lower === upper) return sorted[lower]; + return sorted[lower] + (sorted[upper] - sorted[lower]) * (position - lower); +} + +/** What one penalty scored across every fold. */ +export interface LambdaScore { + lambda: number; + /** Residuals (`actual − predicted`) from validation rows only. */ + residuals: number[]; + meanAbsoluteError: number; +} + +/** + * Score one penalty over the folds. Training folds never contain a row dated + * after their validation fold, by construction of {@link rollingOriginFolds}. + */ +export function scoreLambda( + matrix: FeatureMatrix, + lambda: number, + folds: readonly RollingOriginFold[], +): LambdaScore | null { + const residuals: number[] = []; + for (const fold of folds) { + const fit = fitRidge( + matrix.rows.slice(0, fold.trainEnd), + matrix.targets.slice(0, fold.trainEnd), + lambda, + ); + if (fit === null) continue; + for (let i = fold.validationStart; i < fold.validationEnd; i++) { + const predicted = predictWithFit(fit, matrix.rows[i]); + if (predicted === null) continue; + residuals.push(matrix.targets[i] - predicted); + } + } + if (residuals.length === 0) return null; + return { + lambda, + residuals, + meanAbsoluteError: + residuals.reduce((sum, r) => sum + Math.abs(r), 0) / residuals.length, + }; +} + +/** A model, ready to score a day, with the band its own errors earned it. */ +export interface PrognosisModel { + fit: RidgeFit; + features: string[]; + standardisation: FeatureMatrix["standardisation"]; + /** Rows the final fit used. This is the `n` every surface must show. */ + n: number; + /** Offsets from a prediction to the edges of the band, in target units. */ + bandLow: number; + bandHigh: number; + /** Mean absolute error over the validation folds, for the record. */ + meanAbsoluteError: number; + modelVersion: string; +} + +/** + * Validate, choose λ, refit on the whole window, and take the band from the + * validation residuals. + * + * `null` is a refusal and is a normal outcome: too few complete days to split, + * a singular system, or a fold set that produced no residual at all. The + * caller writes no forecast rather than a forecast nobody measured. + */ +export function fitPrognosisModel( + matrix: FeatureMatrix, + lambdaGrid: readonly number[] = LAMBDA_GRID, +): PrognosisModel | null { + const folds = rollingOriginFolds(matrix.rows.length); + if (folds.length === 0) return null; + + let best: LambdaScore | null = null; + for (const lambda of lambdaGrid) { + const score = scoreLambda(matrix, lambda, folds); + if (score === null) continue; + if (best === null || score.meanAbsoluteError < best.meanAbsoluteError) { + best = score; + } + } + if (best === null) return null; + + const fit = fitRidge(matrix.rows, matrix.targets, best.lambda); + if (fit === null) return null; + + return { + fit, + features: matrix.features, + standardisation: matrix.standardisation, + n: matrix.rows.length, + bandLow: quantile(best.residuals, BAND_LOWER_QUANTILE), + bandHigh: quantile(best.residuals, BAND_UPPER_QUANTILE), + meanAbsoluteError: best.meanAbsoluteError, + modelVersion: MODEL_VERSION, + }; +} diff --git a/src/lib/mood/linked-context.ts b/src/lib/mood/linked-context.ts index f7c858c41..318dcc6cd 100644 --- a/src/lib/mood/linked-context.ts +++ b/src/lib/mood/linked-context.ts @@ -99,6 +99,20 @@ export interface LinkedBody { episodeId: string | null; } +/** + * The five linked figures of one day as plain numbers, for callers that need + * a window of days rather than one day's blocks — the prognosis fit is the + * only one so far. `null` is absence, and a module switched off makes its own + * figures absent rather than zero. + */ +export interface LinkedDayFigures { + sleepAsleep: number | null; + steps: number | null; + activeEnergy: number | null; + restingHeartRate: number | null; + heartRateVariability: number | null; +} + export interface LinkedDayContext { /** The entry's local day, resolved through the entry's own `tz`. */ day: string; @@ -396,6 +410,17 @@ export async function resolveLinkedDayContext( } : MODULE_OFF; + const bodyBlock = await resolveBody(userId, day, enabled); + + return { day, sleep, activity, vitals, body: bodyBlock }; +} + +/** The illness day-log block, split out so the day loop above reads straight. */ +async function resolveBody( + userId: string, + day: string, + enabled: ReadonlyMap, +): Promise> { let body: LinkedBlock = MODULE_OFF; if (enabled.get("illness") ?? true) { // Read only. The illness module owns symptom capture and its severity @@ -419,6 +444,152 @@ export async function resolveLinkedDayContext( episodeId: dayLog?.episodeId ?? null, }; } + return body; +} - return { day, sleep, activity, vitals, body }; +/** + * The five linked figures of one day, as plain numbers, for a whole window. + * + * The same figures the block above answers with, resolved for many days in + * one query instead of one query per day. It exists because the prognosis fits + * over a year of days and the per-day resolver would have run several hundred + * round trips to build one matrix. + * + * It is the SAME engine, deliberately: the module gates come from + * `moduleForMeasurementType` through `typeAvailable`, the cross-source de-dup + * runs through `canonicalRowsOfDay`, the night reconstruction is the one in + * `sleep-night.ts`, and the day totals are the same `sumOfDay` / `latestOfDay` + * the single-day path uses. A second implementation of "what did this day look + * like" would be a second answer, and the mood sheet and the model would start + * disagreeing about the same Tuesday. + * + * A switched-off module answers `null` for its metrics, not zero and not a + * value the owner declined to surface: the model must not learn from a module + * whose figures the account asked not to see. + */ +export async function resolveLinkedDayFigures( + userId: string, + days: readonly string[], + tz: string, +): Promise> { + const out = new Map(); + if (days.length === 0) return out; + + const sorted = [...days].sort(); + const { from } = localDayWindow(sorted[0]); + const { to } = localDayWindow(sorted[sorted.length - 1]); + + const gateStates = await Promise.all( + LINKED_MODULE_KEYS.map( + async (key) => [key, await isModuleEnabled(userId, key)] as const, + ), + ); + const enabled = new Map(gateStates); + const priority = await loadUserSourcePriority(userId); + + const measurements: LinkedRow[] = await prisma.measurement.findMany({ + where: { + userId, + deletedAt: null, + type: { in: [...LINKED_MEASUREMENT_TYPES] }, + measuredAt: { gte: from, lte: to }, + }, + select: { + type: true, + value: true, + measuredAt: true, + sleepStage: true, + source: true, + deviceType: true, + }, + // The canonical-source picker keeps insertion order inside a bucket, so a + // stable read order is what makes its device-type tie-break reproducible + // rather than whatever the planner returned this time. + orderBy: [{ measuredAt: "asc" }, { id: "asc" }], + }); + + // One pass to bucket by local day, rather than filtering the whole window + // once per day: over a year of history that difference is the whole cost of + // the fit's data preparation. + const byDay = new Map(); + for (const row of measurements) { + const key = moodDateKey(row.measuredAt, tz); + const bucket = byDay.get(key); + if (bucket) bucket.push(row); + else byDay.set(key, [row]); + } + + const sleepAvailable = typeAvailable("SLEEP_DURATION", enabled); + const activityAvailable = typeAvailable("ACTIVITY_STEPS", enabled); + const vitalsAvailable = typeAvailable("RESTING_HEART_RATE", enabled); + + const nightByDay = new Map(); + const inBedByDay = new Map(); + if (sleepAvailable) { + const stageRows = measurements.filter((r) => r.type === "SLEEP_DURATION"); + if (stageRows.length > 0) { + for (const night of reconstructSleepNights( + stageRows as unknown as SleepStageRow[], + tz, + priority, + )) { + nightByDay.set(night.night, night.asleepMinutes ?? null); + inBedByDay.set(night.night, night.inBedMinutes ?? null); + } + } + } + + for (const day of days) { + const dayRows = byDay.get(day) ?? []; + out.set(day, { + sleepAsleep: sleepAvailable ? (nightByDay.get(day) ?? null) : null, + steps: activityAvailable + ? valueOf( + sumOfDay(dayRows, "ACTIVITY_STEPS", "steps", day, tz, priority), + ) + : null, + activeEnergy: activityAvailable + ? valueOf( + sumOfDay( + dayRows, + "ACTIVE_ENERGY_BURNED", + "kcal", + day, + tz, + priority, + ), + ) + : null, + restingHeartRate: vitalsAvailable + ? valueOf( + latestOfDay( + dayRows, + "RESTING_HEART_RATE", + "bpm", + day, + tz, + priority, + ), + ) + : null, + heartRateVariability: vitalsAvailable + ? valueOf( + latestOfDay( + dayRows, + "HEART_RATE_VARIABILITY", + "ms", + day, + tz, + priority, + ), + ) + : null, + }); + } + return out; +} + +/** A figure as a plain number, with absence staying absence. */ +function valueOf(figure: LinkedFigure): number | null { + return figure.present ? figure.value : null; }